IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 9627


Ignore:
Timestamp:
Oct 18, 2006, 8:44:12 AM (20 years ago)
Author:
eugene
Message:

select refstars by magnitude, fix various memory leaks, better error tests, fix luminosity function test

Location:
trunk/psastro/src
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/Makefile.am

    r9587 r9627  
    1313libpsastro_la_SOURCES = \
    1414psastroArguments.c          \
     15psastroCleanup.c            \
    1516psastroParseCamera.c        \
    1617psastroDataLoad.c           \

  • trunk/psastro/src/psastro.c

    r9587 r9627  
    22
    33static void usage (void) {
    4     psErrorStackPrint(stderr, "USAGE: psastro [-file image(s)] [-list imagelist] (output)");
     4    psErrorStackPrint(stderr, "USAGE: psastro [-file image(s)] [-list imagelist] (output)\n");
    55    exit (2);
    66}
    … …  
    2020    // load identify the data sources
    2121    if (!psastroParseCamera (config)) {
    22         psErrorStackPrint(stderr, "error setting up the camera");
     22        psErrorStackPrint(stderr, "error setting up the camera\n");
    2323        exit (1);
    2424    }
    … …  
    2727    // select subset of stars for astrometry
    2828    if (!psastroDataLoad (config)) {
    29         psErrorStackPrint(stderr, "error loading input data");
     29        psErrorStackPrint(stderr, "error loading input data\n");
    3030        exit (1);
    3131    }
    … …  
    3434    // apply the initial guess
    3535    if (!psastroAstromGuess (config)) {
    36         psErrorStackPrint(stderr, "failed to determine initial astrometry guess");
     36        psErrorStackPrint(stderr, "failed to determine initial astrometry guess\n");
    3737        exit (1);
    3838    }
    … …  
    4141    psArray *refs = psastroLoadRefstars (config);
    4242    if (!refs) {
    43         psErrorStackPrint(stderr, "failed to load reference data");
     43        psErrorStackPrint(stderr, "failed to load reference data\n");
    4444        exit (1);
    4545    }
    4646
    4747    // choose reference stars corresponding to the selected chips
    48     if (psastroChooseRefstars (config, refs)) {
    49         psErrorStackPrint(stderr, "failed to select reference data for chips");
     48    if (!psastroChooseRefstars (config, refs)) {
     49        psErrorStackPrint(stderr, "failed to select reference data for chips\n");
    5050        exit (1);
    5151    }
    … …  
    5454    if (mosastro == NULL) {
    5555        if (!psastroChipAstrom (config, refs)) {
    56             psErrorStackPrint(stderr, "failed to perform single chip astrometry");
     56            psErrorStackPrint(stderr, "failed to perform single chip astrometry\n");
    5757            exit (1);
    5858        }
    5959    } else {
    6060        if (!psastroMosaicAstrom (config, refs)) {
    61             psErrorStackPrint(stderr, "failed to perform mosaic camera astrometry");
     61            psErrorStackPrint(stderr, "failed to perform mosaic camera astrometry\n");
    6262            exit (1);
    6363        }
    … …  
    6969    // write out the results
    7070    if (!psastroDataSave (config)) {
    71         psErrorStackPrint(stderr, "failed to write out data");
     71        psErrorStackPrint(stderr, "failed to write out data\n");
    7272        exit (1);
    7373    }
    … …  
    7575    psLogMsg ("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete"));
    7676
    77     psFree (refs);
    78     psFree (config);
    79 
    80     psTimerStop ();
    81     psMemCheckCorruption (true);
    82     pmModelGroupCleanup ();
    83     psTimeFinalize ();
    84     pmConceptsDone ();
    85     pmConfigDone ();
    86     // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro");
    87     fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro");
    88 
    89     exit (0);
     77    psastroCleanup (config, refs);
     78    exit (EXIT_SUCCESS);
    9079}

  • trunk/psastro/src/psastro.h

    r9587 r9627  
    3131
    3232pmConfig         *psastroArguments (int argc, char **argv);
     33void              psastroCleanup (pmConfig *config, psArray *refs);
    3334bool              psastroParseCamera (pmConfig *config);
    3435bool              psastroDataLoad (pmConfig *config);
    … …  
    4445psPlaneDistort   *psPlaneDistortIdentity ();
    4546
    46 psArray          *psastroLoadReferences (pmConfig *config);
     47psArray          *psastroLoadRefstars (pmConfig *config);
    4748bool              psastroChipAstrom (pmConfig *config, psArray *refs);
    4849bool              psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refset, psArray *rawset, psMetadata *recipe, psMetadata *updates);
     50bool              psastroChooseRefstars (pmConfig *config, psArray *refs);
     51bool              psastroRefstarSubset (pmReadout *readout);
     52pmLumFunc        *psastroLuminosityFunction (psArray *stars);
    4953
    5054// utility functions:
    … …  
    5660psArray *psastroMosaicGetGrads (pmFPA *fpa, psMetadata *recipe);
    5761bool psastroMosaicAstrom (pmConfig *config, psArray *refs);
    58 bool psastroMosaicGetRefstars (pmConfig *config, psArray *refs);
    5962bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe);
    6063bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe);

  • trunk/psastro/src/psastroArguments.c

    r9574 r9627  
    2727    // load config data from default locations
    2828    pmConfig *config = pmConfigRead(&argc, argv);
     29    if (config == NULL) {
     30        psError(PSASTRO_ERR_CONFIG, false, "Can't read site configuration");
     31        return NULL;
     32    }
    2933
    3034    // save these recipe options until we have loaded the options

  • trunk/psastro/src/psastroAstromGuess.c

    r9574 r9627  
    8080                    p_psDeproject (raw->sky, raw->TP, fpa->projection);
    8181
    82                     if ((i < 10) && (psTraceGetLevel("psastro.guess") > 5)) {
     82                    if ((i < 10) && (psTraceGetLevel("psastro") > 5)) {
    8383                        fprintf (stderr, "up: %f,%f -> %f,%f -> %f,%f -> %f,%f\n",
    8484                                 raw->chip->x, raw->chip->y,
    … …  
    100100                                 tp->x, tp->y,
    101101                                 raw->sky->r, raw->sky->d);
     102
     103                        psFree (fp);
     104                        psFree (tp);
     105                        psFree (ch);
    102106                    }
    103107

  • trunk/psastro/src/psastroChipAstrom.c

    r9574 r9627  
    5454                pmAstromWriteWCS (chip->toFPA, fpa->toSky, updates, plateScale);
    5555                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_DATA_METADATA, "psastro header stats", updates);
     56                psFree (updates);
    5657            }
    5758        }

  • trunk/psastro/src/psastroChooseRefstars.c

    r9587 r9627  
    7979
    8080                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
     81                psFree (refstars);
    8182
    8283                psastroRefstarSubset (readout);
    83 
    84                 psFree (refstars);
    8584            }
    8685        }

  • trunk/psastro/src/psastroErrorCodes.dat

    r9575 r9627  
    77ARGUMENTS               Incorrect arguments
    88CONFIG                  Problem in configure files
    9 FITS                    Problem in FITS I/O
    10 FITS_WCS                Error interpreting FITS WCS information
    11 PHOTOM                  Problem in photometry
    12 PSF                     Problem in PSF
    13 APERTURE                Problem with aperture photometry
     9IO                      Problem in FITS I/O
     10WCS                     Error interpreting FITS WCS information
     11DATA                    Problem in data values

  • trunk/psastro/src/psastroLoadRefstars.c

    r9587 r9627  
    2020    float DECmin = DEG_RAD*psMetadataLookupF32(NULL, recipe, "DEC_MIN");
    2121    float DECmax = DEG_RAD*psMetadataLookupF32(NULL, recipe, "DEC_MAX");
     22    float MAGmax = psMetadataLookupF32(NULL, recipe, "MAG_MAX");
    2223
    2324    char *CATDIR = psMetadataLookupStr(NULL, recipe, "DVO.CATDIR");
    … …  
    3839    close (fd);
    3940
     41    psTimerStart ("psastro");
     42
    4043    // use fork to add timeout capability
    4144    // use cfitsio |filename format to avoid tempFile
    4245    # if ELIXIR_MODE
    4346    if (CATDIR) {
    44         sprintf (tempLine, "getstar -D CATDIR %s -D CATFORMAT elixir -D CATMODE mef -region %f %f %f %f -o %s", CATDIR, RAmin, DECmin, RAmax, DECmax, tempFile);
     47        sprintf (tempLine, "getstar -D CATDIR %s -D CATFORMAT elixir -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", CATDIR, MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile);
    4548    } else {
    46         sprintf (tempLine, "getstar -D CATFORMAT elixir -D CATMODE mef -region %f %f %f %f -o %s", RAmin, DECmin, RAmax, DECmax, tempFile);
     49        sprintf (tempLine, "getstar -D CATFORMAT elixir -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile);
    4750    }
    4851    # else
    4952    if (CATDIR) {
    50         sprintf (tempLine, "getstar -D CATDIR %s -D CATFORMAT panstarrs -D CATMODE mef -region %f %f %f %f -o %s", CATDIR, RAmin, DECmin, RAmax, DECmax, tempFile);
     53        sprintf (tempLine, "getstar -D CATDIR %s -D CATFORMAT panstarrs -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", CATDIR, MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile);
    5154    } else {
    52         sprintf (tempLine, "getstar -D CATFORMAT panstarrs -D CATMODE mef -region %f %f %f %f -o %s", RAmin, DECmin, RAmax, DECmax, tempFile);
     55        sprintf (tempLine, "getstar -D CATFORMAT panstarrs -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile);
    5356    }
    5457    # endif
    … …  
    6164    }
    6265
     66    psLogMsg ("psastro", 3, "ran getstar : %f sec\n", psTimerMark ("psastro"));
     67
    6368    // the output from getstar is a file with the Average table
    6469    psFits *fits = psFitsOpen (tempFile, "r");
    … …  
    7075    # endif
    7176
     77    psTimerStart ("psastro");
    7278    psMetadata *header = psFitsReadHeader (NULL, fits);
    7379    psArray *table = psFitsReadTable (fits);
    7480    psFitsClose (fits);
     81
    7582    unlink (tempFile);
     83    psLogMsg ("psastro", 3, "read getstar output table : %f sec\n", psTimerMark ("psastro"));
    7684
    7785    // convert the Average table to the pmAstromObj entries
     86    psTimerStart ("psastro");
    7887    psArray *refs = psArrayAlloc (table->n);
    7988    for (int i = 0; i < table->n; i++) {
    … …  
    98107    psFree (header);
    99108    psFree (table);
     109    psLogMsg ("psastro", 3, "converted table to pmAstromObj : %f sec\n", psTimerMark ("psastro"));
    100110
    101111    psTrace ("psastro", 3, "loaded %ld reference stars from (%10.6f,%10.6f) - (%10.6f,%10.6f)\n",
    102112             refs->n, RAmin, DECmin, RAmax, DECmax);
     113
    103114    return refs;
    104115}

  • trunk/psastro/src/psastroLuminosityFunction.c

    r9587 r9627  
    1515  psMemSetDeallocator(func, (psFreeFunc) pmLumFuncFree);
    1616
    17   func->mMin = 0;
    18   func->mMax = 0;
    19   func->slope = 0;
    20   func->offset = 0;
     17  func->mMin = mMin;
     18  func->mMax = mMax;
     19  func->slope = slope;
     20  func->offset = offset;
    2121
    2222  return func;
    2323}
    2424
    25 pmLumFunc *psastroRefstarSubset (psArray *stars) {
     25pmLumFunc *psastroLuminosityFunction (psArray *stars) {
    2626
    2727  // determine the max and min magnitude for the array of stars
    2828  pmAstromObj *star = stars->data[0];
    29   mMin = star->Mag;
    30   mMax = star->Mag;
     29  double mMin = star->Mag;
     30  double mMax = star->Mag;
    3131  for (int i = 0; i < stars->n; i++) {
    3232    star = stars->data[i];
    … …  
    3535  }
    3636
    37   int nBin = 1 + (nMax - nMin) / dMag;
     37  int nBin = 1 + (mMax - mMin) / dMag;
    3838  if (nBin <= 1) {
    3939    psError(PSASTRO_ERR_DATA, true, "magnitude range of 0.0\n");
    … …  
    4545  // bin[i] = mMin + i*dMag : mMin + (i+1)*dMag
    4646  psVector *nMags = psVectorAlloc (nBin, PS_TYPE_F32);
     47  nMags->n = nBin;
     48  psVectorInit (nMags, 0);
     49
    4750  for (int i = 0; i < stars->n; i++) {
    4851    star = stars->data[i];
    4952    int bin = (star->Mag - mMin) / dMag;
    5053    nMags->data.F32[bin] += 1.0;
     54  }
     55
     56  // find the peak and position
     57  int iPeak = 0;
     58  int nPeak = 0;
     59  for (int i = 0; i < nMags->n; i++) {
     60    if (nMags->data.F32[i] < nPeak) continue;
     61    iPeak = i;
     62    nPeak = nMags->data.F32[i];
    5163  }
    5264
    … …  
    5668  // create 2 vectors represnting the dlogN/dlogS line
    5769  // exclude bins with no stars
     70  // exclude points after the peak with N < 0.8*peak
    5871  int n = 0;
    59   for (int i = 0; i < stars->n; i++) {
     72  for (int i = 0; i < nMags->n; i++) {
    6073    if (nMags->data.F32[i] < 1) continue;
     74    if ((i > iPeak) && (nMags->data.F32[i] < 0.8*nPeak)) continue;
    6175    lnMag->data.F32[n] = log10 (nMags->data.F32[i]);
    6276    Mag->data.F32[n] = mMin + (i + 0.5)*dMag;
    6377    n++;
    6478  }
     79  lnMag->n = n;
     80  Mag->n = n;
     81  psLogMsg ("psastro", 4, "fitting %d points to luminosity function\n", n);
    6582
    6683  psVector *mask = psVectorAlloc (nBin, PS_TYPE_MASK);
     84  mask->n = Mag->n;
     85  psVectorInit (mask, 0);
     86
    6787  psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
    6888  psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    … …  
    7090  stats->clipIter = 3;
    7191
    72   line = psVectorClipFitPolynomial1D(line, stats, mask, 1, lnMag, NULL, Mag);
     92  line = psVectorClipFitPolynomial1D(line, stats, mask, 0xff, lnMag, NULL, Mag);
    7393
    7494  // find min and max unmasked magnitudes
    75   mMinValid = NAN;
    76   mMaxValid = NAN;
     95  double mMinValid = NAN;
     96  double mMaxValid = NAN;
    7797  for (int i = 0; i < Mag->n; i++) {
    7898    if (mask->data.U8[i]) continue;
    … …  
    89109
    90110  pmLumFunc *lumFunc = pmLumFuncAlloc (mMinValid, mMaxValid, line->coeff[0], line->coeff[1]);
     111
     112  psFree (lnMag);
     113  psFree (nMags);
     114  psFree (Mag);
     115  psFree (mask);
     116  psFree (line);
     117  psFree (stats);
     118
    91119  return (lumFunc);
    92120}

  • trunk/psastro/src/psastroRefstarSubset.c

    r9587 r9627  
    55  // select the raw objects for this readout
    66  psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    7   if (rawstars == NULL) { continue; }
     7  if (rawstars == NULL)  {
     8    psError(PSASTRO_ERR_DATA, false, "missing rawstars in psastroRefstarSubset\n");
     9    return false;
     10  }
    811
    912  // select the raw objects for this readout
    1013  psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    11   if (refstars == NULL) { continue; }
     14  if (refstars == NULL)  {
     15    psError(PSASTRO_ERR_DATA, false, "missing refstars in psastroRefstarSubset\n");
     16    return false;
     17  }
    1218
    1319  // calculate luminosity functions for rawstars and refstars
    … …  
    1521  // is needed...
    1622  pmLumFunc *rawfunc = psastroLuminosityFunction (rawstars);
     23  if (rawfunc == NULL) {
     24    psError(PSASTRO_ERR_DATA, false, "error measuring rawstar luminosity function\n");
     25    return NULL;
     26  }
    1727  pmLumFunc *reffunc = psastroLuminosityFunction (refstars);
     28  if (reffunc == NULL) {
     29    psError(PSASTRO_ERR_DATA, false, "error measuring refstar luminosity function\n");
     30    return NULL;
     31  }
     32 
     33  // XXXX test
     34  // psFree (rawfunc);
     35  // psFree (reffunc);
     36  // return true;
    1837
    1938  // what is the offset between the two lines at the average magnitude?
    20   mRef = 0.5*(reffunc->mMin + reffunc->mMax);
    21   logRho = mRef * reffunc->slope + reffunc->offset;
    22   mRaw = (logRho - rawfunc->offset) / rawfunc->slope;
     39  double mRef = 0.5*(reffunc->mMin + reffunc->mMax);
     40  double logRho = mRef * reffunc->slope + reffunc->offset;
     41  double mRaw = (logRho - rawfunc->offset) / rawfunc->slope;
    2342
    2443  psLogMsg ("psastro", 4, "mRef: %f, logRho: %f, mRaw: %f\n", mRef, logRho, mRaw);
    2544
    26   mRefMax = rawfunc->mMax - mRaw + mRef;
     45  double mRefMax = rawfunc->mMax - mRaw + mRef;
    2746  psLogMsg ("psastro", 4, "clipping stars fainter than %f\n", mRefMax);
    2847
    … …  
    3150    pmAstromObj *ref = refstars->data[i];
    3251    if (ref->Mag > mRefMax) continue;
    33     psArrayAdd (subset, ref);
     52    psArrayAdd (subset, 100, ref);
    3453  }
    3554
    36   psLogMsg ("psastro", 4, "keeping %d of %d reference stars\n", subset->n, refstars->n);
     55  psLogMsg ("psastro", 4, "keeping %ld of %ld reference stars\n", subset->n, refstars->n);
    3756
    3857  psMetadataRemoveKey (readout->analysis, "PSASTRO.REFSTARS");
    3958  psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", subset);
    4059
     60  psFree (rawfunc);
     61  psFree (reffunc);
     62  psFree (subset);
     63
    4164  return true;
    4265}
Note: See TracChangeset for help on using the changeset viewer.