IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 26, 2009, 1:59:32 PM (17 years ago)
Author:
beaumont
Message:

merged with trunk

Location:
branches/cnb_branches/cnb_branch_20090301
Files:
1 deleted
26 edited
6 copied

Legend:

Unmodified
Added
Removed

  • branches/cnb_branches/cnb_branch_20090301

  • branches/cnb_branches/cnb_branch_20090301/psastro

  • branches/cnb_branches/cnb_branch_20090301/psastro/src

    • Property svn:ignore

      •  

        old new  
        1717psastroModelFit
        1818psastroExtract
         19psastroVersionDefinitions.h

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/Makefile.am

    r23594 r24244  
    11lib_LTLIBRARIES = libpsastro.la
    22
    3 libpsastro_la_CFLAGS = $(PSASTRO_CFLAGS) $(PPSTATS_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) -DPSASTRO_VERSION=$(SVN_VERSION) -DPSASTRO_BRANCH=$(SVN_BRANCH) -DPSASTRO_SOURCE=$(SVN_SOURCE)
    4 libpsastro_la_LDFLAGS = $(PSASTRO_LIBS) $(PPSTATS_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     3if HAVE_SVNVERSION
     4PSASTRO_VERSION=`$(SVNVERSION) ..`
     5else
     6PSASTRO_VERSION="UNKNOWN"
     7endif
     8
     9if HAVE_SVN
     10PSASTRO_BRANCH=`$(SVN) info .. | $(SED) -n -e '/URL:/ h' -e '/Repository Root:/ { x; H; x; s|Repository Root: \(.*\)\nURL: \1\(.*\)|\2| ; s|^/|| ; s|/[a-zA-Z]*/src.*|| ; p }'`
     11PSASTRO_SOURCE=`$(SVN) info | $(SED) -n -e 's/Repository UUID: // p'`
     12else
     13PSASTRO_BRANCH="UNKNOWN"
     14PSASTRO_SOURCE="UNKNOWN"
     15endif
    516
    617# Force recompilation of psastroVersion.c, since it gets the version information
    7 # can we do this with dependency info?
    8 # psastroVersion.c: FORCE
    9 # touch psastroVersion.c
    10 # FORCE: ;
     18psastroVersion.c: psastroVersionDefinitions.h
     19psastroVersionDefinitions.h: psastroVersionDefinitions.h.in FORCE
     20        -$(RM) psastroVersionDefinitions.h
     21        $(SED) -e "s|@PSASTRO_VERSION@|\"$(PSASTRO_VERSION)\"|" -e "s|@PSASTRO_BRANCH@|\"$(PSASTRO_BRANCH)\"|" -e "s|@PSASTRO_SOURCE@|\"$(PSASTRO_SOURCE)\"|" psastroVersionDefinitions.h.in > psastroVersionDefinitions.h
     22FORCE: ;
     23
     24libpsastro_la_CFLAGS = $(PSASTRO_CFLAGS) $(PPSTATS_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     25libpsastro_la_LDFLAGS = $(PSASTRO_LIBS) $(PPSTATS_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    1126
    1227bin_PROGRAMS = psastro psastroExtract psastroModel psastroModelFit gpcModel
    … …  
    5065        psastroExtractStars.c      \
    5166        psastroExtractGhosts.c      \
    52         psastroExtractFindChip.c      \
    5367        psastroCleanup.c
    5468
    … …  
    7387
    7488libpsastro_la_SOURCES = \
     89        psastroMaskUpdates.c        \
     90        psastroMaskUtils.c          \
     91        psastroLoadGhosts.c         \
     92        psastroGhostUtils.c         \
    7593        psastroErrorCodes.c         \
    7694        psastroVersion.c            \
    7795        psastroDefineFiles.c        \
    7896        psastroAnalysis.c           \
    79         psastroMaskUpdates.c        \
    8097        psastroAstromGuess.c        \
    8198        psastroLoadRefstars.c       \
    … …  
    102119        psastroMosaicSetAstrom.c    \
    103120        psastroMosaicSetMatch.c     \
     121        psastroFindChip.c           \
    104122        psastroZeroPoint.c          \
    105123        psastroDemoDump.c           \
    … …  
    122140
    123141### Error codes.
    124 BUILT_SOURCES = psastroErrorCodes.h psastroErrorCodes.c
     142BUILT_SOURCES = psastroErrorCodes.h psastroErrorCodes.c psastroVersionDefinitions.h
    125143CLEANFILES = psastroErrorCodes.h psastroErrorCodes.c
    126144EXTRA_DIST = psastroErrorCodes.dat psastroErrorCodes.h.in psastroErrorCodes.c.in

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastro.c

    r23352 r24244  
    1313# include "psastroStandAlone.h"
    1414
    15 static void usage (void) {
    16     fprintf (stderr, "USAGE: psastro [-file image(s)] [-list imagelist] (output)\n");
    17     exit (2);
     15static void usage(void) {
     16    fprintf(stderr, "USAGE: psastro [-file image(s)] [-list imagelist] (output)\n");
     17    exit(PS_EXIT_CONFIG_ERROR);
    1818}
    1919
    20 int main (int argc, char **argv) {
    21 
    22     pmConfig *config = NULL;
    23 
     20int main (int argc, char **argv)
     21{
    2422    psTimerStart ("complete");
    2523
    … …  
    2826    // model inits are needed in pmSourceIO
    2927    // models defined in psphot/src/models are not available in psastro
    30     pmModelClassInit ();
     28    pmModelClassInit();
    3129
    3230    // load configuration information
    33     config = psastroArguments (argc, argv);
    34     if (!config) usage ();
     31    pmConfig *config = config = psastroArguments(argc, argv);
     32    if (!config) {
     33        usage();
     34    }
    3535
    3636    psastroVersionPrint();
    … …  
    3939    if (!psastroParseCamera (config)) {
    4040        psErrorStackPrint(stderr, "error setting up the camera\n");
    41         exit (1);
     41        psFree(config);
     42        exit(PS_EXIT_CONFIG_ERROR);
    4243    }
    4344
    … …  
    4647    if (!psastroDataLoad (config)) {
    4748        psErrorStackPrint(stderr, "error loading input data\n");
    48         exit (1);
     49        psFree(config);
     50        exit(PS_EXIT_DATA_ERROR);
    4951    }
    5052
     53    psMetadata *stats = psMetadataAlloc(); // Statistics, for output
     54    psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0);
     55
    5156    // run the full astrometry analysis (chip and/or mosaic)
    52     if (!psastroAnalysis (config)) {
     57    if (!psastroAnalysis(config, stats)) {
    5358        psErrorStackPrint(stderr, "failure in psastro analysis\n");
    54         exit (1);
     59        psFree(config);
     60        psFree(stats);
     61        exit(PS_EXIT_SYS_ERROR);
    5562    }
    5663
    5764    // write out the results
    58     if (!psastroDataSave (config)) {
     65    if (!psastroDataSave(config, stats)) {
    5966        psErrorStackPrint(stderr, "failed to write out data\n");
    60         exit (1);
     67        psFree(config);
     68        psFree(stats);
     69        exit(PS_EXIT_DATA_ERROR);
    6170    }
    6271
    63     psLogMsg ("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete"));
     72    psFree (stats);
     73    psLogMsg("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete"));
    6474
    65     psastroCleanup (config);
    66     exit (EXIT_SUCCESS);
     75    psastroCleanup(config);
     76    exit(PS_EXIT_SUCCESS);
    6777}

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastro.h

    r23594 r24244  
    4646#endif
    4747
    48 bool              psastroDataSave (pmConfig *config);
     48/**
     49 * this structure defines the parameters to describe a ghost
     50 */
     51typedef struct {
     52    psPlane *srcFP;                     ///< location in FPA coords of the source star
     53    psPlane *FP;                        ///< location in FPA coords of the ghost center
     54    psPlane *chip;                      ///< location in chip coords of the ghost center
     55    double Mag;                         ///< instrumental magnitude of source star
     56    psEllipseAxes inner;                ///< inner elliptical annulus boundary
     57    psEllipseAxes outer;                ///< outer elliptical annulus boundary
     58} psastroGhost;
     59
     60psastroGhost     *psastroGhostAlloc (void);
     61bool              psastroLoadGhosts (pmConfig *config);
     62
     63bool              psastroDataSave (pmConfig *config, psMetadata *stats);
    4964bool              psastroDefineFiles (pmConfig *config, pmFPAfile *input);
    5065bool              psastroDefineFile (pmConfig *config, pmFPA *input, char *filerule, char *argname, pmFPAfileType fileType, pmDetrendType detrendType);
    51 bool              psastroAnalysis (pmConfig *config);
     66bool              psastroAnalysis (pmConfig *config, psMetadata *stats);
    5267
    5368bool              psastroConvertFPA (pmFPA *fpa, psMetadata *recipe);
    … …  
    125140bool              psastroAstromGuessSetChip (pmFPA *fpa, pmChip *chip, const pmFPAview *view, double pixelScale, bool bilevelAstrometry);
    126141bool              psastroAstromGuessSetFPA (pmFPA *fpa, bool *bilevelAstrometry);
    127 bool              psastroMetadataStats (pmConfig *config);
     142bool              psastroMetadataStats (pmConfig *config, psMetadata *stats);
    128143
    129144bool              psastroZeroPoint (pmConfig *config);
    … …  
    131146bool              psastroZeroPointFromRecipe (float *zeropt, float *exptime, pmFPA *fpa, psMetadata *recipe);
    132147
     148// masking functions
     149pmCell           *pmCellInChip (pmChip *chip, float x, float y);
     150bool              pmCellCoordsForChip (float *xCell, float *yCell, pmCell *cell, float xChip, float yChip);
     151bool              pmChipCoordsForCell (float *xChip, float *yChip, pmCell *cell, float xCell, float yCell);
     152                 
     153bool              psastroMaskCircle (psImage *mask, psImageMaskType value, float x0, float y0, float dX, float dY);
     154bool              psastroMaskEllipse (psImage *mask, psImageMaskType value, float x0, float y0, psEllipseAxes axes);
     155bool              psastroMaskEllipticalAnnulus (psImage *mask, psImageMaskType value, float x0, float y0, psEllipseAxes eInner, psEllipseAxes eOuter);
     156                 
     157bool              psastroMaskBox (psImage *mask, psImageMaskType value, float x0, float y0, float dL, float dW, float theta);
     158void              psastroMaskLine (psImage *mask, psImageMaskType value, double x1, double y1, double x2, double y2, int dW);
     159void              psastroMaskLineBresen (psImage *mask, psImageMaskType value, int X1, int Y1, int X2, int Y2, int dW, int swapcoords);
     160void              psastroMaskRectangle (psImage *mask, psImageMaskType value, int x0, int y0, int x1, int y1);
     161
     162pmChip           *psastroFindChip (double *xChip, double *yChip, pmFPA *fpa, double xFPA, double yFPA);
     163bool              psastroChipBounds (pmFPA *fpa);
     164
    133165// psastroExtract functions
    134166bool              psastroExtractAnalysis (pmConfig *config);
    135167bool              psastroExtractStars (pmConfig *config);
    136168bool              psastroExtractGhosts (pmConfig *config);
    137 pmChip           *psastroExtractFindChip (float *xChip, float *yChip, pmFPA *fpa, float xFPA, float yFPA);
    138 bool              psastroExtractChipBounds (pmFPA *fpa);
    139169
    140 psImage          *psastroExtractStar (psImage *input, float x, float y, float dX, float dY);
     170psImage          *psastroExtractStar (psImage *input, double x, double y, double dX, double dY);
    141171bool              psastroExtractParseCamera (pmConfig *config);
    142172bool              psastroExtractDataLoad (pmConfig *config);
    143173pmConfig         *psastroExtractArguments (int argc, char **argv);
     174bool              psastroExtractFreeChipBounds(void);
    144175
    145176///@}

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroAnalysis.c

    r23594 r24244  
    11/** @file psastroAnalysis.c
    22 *
    3  *  @brief 
     3 *  @brief
    44 *
    55 *  @ingroup libpsastro
    … …  
    1313# include "psastroInternal.h"
    1414
    15 bool psastroAnalysis (pmConfig *config) {
     15/// Turn save on/off for a file
     16static void fileSave(pmConfig *config,  // Configuration
     17                     const char *name,  // Name of file
     18                     bool save          // Save file?
     19    )
     20{
     21    pmFPAfile *file = pmFPAfileSelectSingle(config->files, name, 0); // File of interest
     22    if (!file) {
     23        psErrorClear();
     24        return;
     25    }
     26    file->save = save;
     27    return;
     28}
     29
     30
     31bool psastroAnalysis (pmConfig *config, psMetadata *stats) {
    1632
    1733    bool status;
    … …  
    4056    }
    4157    if (nStars == 0) {
    42         psLogMsg ("psastro", 2, "skipping astrometry analysis : no stars\n");
    43         return false;
     58        // This is likely a data quality issue
     59        psWarning("No stars for astrometry analysis --- suspect bad data quality");
     60        if (stats && psMetadataLookupS32(NULL, stats, "QUALITY") == 0) {
     61            psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
     62                             "No stars for astrometry", PSASTRO_ERR_DATA);
     63        }
     64        fileSave(config, "PSASTRO.OUTPUT", false);
     65        fileSave(config, "PSASTRO.OUTPUT.MASK", false);
     66        fileSave(config, "PSASTRO.OUT.REFSTARS", false);
     67        psErrorClear();
     68        return true;
    4469    }
    4570

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroAstromGuess.c

    r21422 r24244  
    382382    psStats *statsQ = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    383383
    384     psVectorStats (statsP, cornerPd, NULL, cornerMK, 1);
    385     psVectorStats (statsQ, cornerQd, NULL, cornerMK, 1);
     384    if (!psVectorStats (statsP, cornerPd, NULL, cornerMK, 1)) {
     385        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     386        return false;
     387    }
     388    if (!psVectorStats (statsQ, cornerQd, NULL, cornerMK, 1)) {
     389        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     390        return false;
     391    }
    386392
    387393    float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]);

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroChipAstrom.c

    r21409 r24244  
    130130    }
    131131
     132    if (numGoodChips == 0) {
     133        psError(PSASTRO_ERR_UNKNOWN, false, "Failed to fit any chips");
     134        return false;
     135    }
     136
    132137    if (!psastroFixChips (config, recipe)) {
    133138        psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips");

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroDataSave.c

    r23352 r24244  
    2222 * this loop saves the photometry/astrometry data files
    2323 */
    24 bool psastroDataSave (pmConfig *config) {
     24bool psastroDataSave (pmConfig *config, psMetadata *stats) {
    2525
    2626    pmChip *chip;
    … …  
    8484
    8585    // Write out summary statistics
    86     if (!psastroMetadataStats (config)) ESCAPE;
     86    if (!psastroMetadataStats (config, stats)) ESCAPE;
     87
     88    bool status;
     89    psString dump_file = psMetadataLookupStr(&status, config->arguments, "DUMP_CONFIG");
     90    if (dump_file) {
     91        pmConfigCamerasCull(config, NULL);
     92        pmConfigRecipesCull(config, "PPIMAGE,PPSTATS,PSPHOT,MASKS,PSASTRO");
     93
     94        pmConfigDump(config, dump_file);
     95    }
    8796
    8897    // activate all files except PSASTRO.OUTPUT

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroErrorCodes.dat

    r11269 r24244  
    22# This file is used to generate pmErrorClasses.h
    33#
    4 BASE = 300              First value we use; lower values belong to psLib
    5 UNKNOWN                 Unknown PM error code
     4BASE = 4000             First value we use; lower values belong to psLib
     5UNKNOWN                 Unknown psastro error code
    66NOT_IMPLEMENTED         Desired feature is not yet implemented
    77ARGUMENTS               Incorrect arguments

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroErrorCodes.h.in

    r21409 r24244  
    2121 */
    2222typedef enum {
    23     PSASTRO_ERR_BASE = 600,
     23    PSASTRO_ERR_BASE = 4000,
    2424    PSASTRO_ERR_${ErrorCode},
    2525    PSASTRO_ERR_NERROR

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroExtractGhosts.c

    r23594 r24244  
    8989
    9090                    // XXX can make a more clever model...
    91                     float xFPA = -1*ref->FP->x;
    92                     float yFPA = -1*ref->FP->y;
     91                    double xFPA = -1*ref->FP->x;
     92                    double yFPA = -1*ref->FP->y;
    9393
    94                     float xChip, yChip;
    95                     pmChip *ghostChip = psastroExtractFindChip (&xChip, &yChip, fpa, xFPA, yFPA);
     94                    double xChip, yChip;
     95                    pmChip *ghostChip = psastroFindChip (&xChip, &yChip, fpa, xFPA, yFPA);
    9696                    if (!ghostChip) continue;
    9797                    if (!ghostChip->cells) continue;
    … …  
    107107                    if (ghosts == NULL) {
    108108                        ghosts = psArrayAllocEmpty (100);
    109                         psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.GHOSTS", PS_DATA_ARRAY, "astrometry matches", ghosts);
     109                        if (!psMetadataAdd (ghostReadout->analysis, PS_LIST_TAIL, "PSASTRO.GHOSTS", PS_DATA_ARRAY, "astrometry matches", ghosts)) {
     110                          psError(PSASTRO_ERR_CONFIG, false, "failure to add ghosts to readout");
     111                          return false;
     112                        }
    110113                        psFree (ghosts);
    111114                    }
    … …  
    127130    }
    128131
     132    psastroExtractFreeChipBounds();
     133
    129134    psFree (view);
    130135    return true;

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroExtractStar.c

    r23594 r24244  
    1313# include "psastroInternal.h"
    1414
    15 psImage *psastroExtractStar (psImage *input, float x, float y, float dX, float dY) {
     15psImage *psastroExtractStar (psImage *input, double x, double y, double dX, double dY) {
    1616
    1717    // skip if no pixels are on the image

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroExtractStars.c

    r23594 r24244  
    7676
    7777    int nExt = 0;
     78    int nExtGhost = 0;
    7879
    7980    // open the output file handle: we are just saving a series of extensions to this file
    … …  
    157158                    psMetadataAddF32 (header, PS_LIST_TAIL, "FPA_Y",   0, "focal-plane coordinate", ref->FP->y);
    158159                    psMetadataAddF32 (header, PS_LIST_TAIL, "MAG",     0, "magnitude",              ref->Mag);
     160                    psMetadataAddF32 (header, PS_LIST_TAIL, "INST_MAG",0, "instrumental magnitude", ref->Mag - MagOffset);
    159161
    160162                    snprintf (extname, 80, "extname.%05d", nExt);
    … …  
    175177
    176178                    pmReadout *inReadout = pmFPAviewThisReadout (view, input->fpa);
    177                     psImage *subraster = psastroExtractStar (inReadout->image, ghost->chip->x, ghost->chip->y, 200.0, 200.0);
     179                    psImage *subraster = psastroExtractStar (inReadout->image, ghost->chip->x, ghost->chip->y, 400.0, 400.0);
    178180                    if (!subraster) continue;
    179181                   
    … …  
    187189                    psMetadataAddF32 (header, PS_LIST_TAIL, "FPA_Y0",  0, "star focal-plane coordinate", ghost->TP->y); // note : over-loaded value
    188190                    psMetadataAddF32 (header, PS_LIST_TAIL, "MAG",     0, "magnitude",              ghost->Mag);
    189 
    190                     snprintf (extname, 80, "extname.%05d", nExt);
    191                     psFitsWriteImage (outStars, header, subraster, 0, extname);
    192                     nExt ++;
     191                    psMetadataAddF32 (header, PS_LIST_TAIL, "INST_MAG",0, "instrumental magnitude", ghost->Mag - MagOffset);
     192
     193                    snprintf (extname, 80, "extname.%05d", nExtGhost);
     194                    psFitsWriteImage (outGhosts, header, subraster, 0, extname);
     195                    nExtGhost ++;
    193196                   
    194197                    psFree (header);
    … …  
    204207
    205208    psFitsClose (outStars);
     209    psFitsClose (outGhosts);
    206210    psFree (view);
    207211    return true;

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroFixChipsTest.c

    r21409 r24244  
    140140    psStats *stats;
    141141    stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    142     psVectorStats (stats, dX, NULL, NULL, 0);
     142    if (!psVectorStats (stats, dX, NULL, NULL, 0)) {
     143        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     144        return false;
     145    }
    143146    Xo = stats->robustMedian;
    144147    fprintf (stderr, "offset x: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
    … …  
    146149
    147150    stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    148     psVectorStats (stats, dY, NULL, NULL, 0);
     151    if (!psVectorStats (stats, dY, NULL, NULL, 0)) {
     152        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     153        return false;
     154    }
    149155    Yo = stats->robustMedian;
    150156    fprintf (stderr, "offset y: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
    … …  
    168174    }
    169175    stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    170     psVectorStats (stats, dT, NULL, NULL, 0);
     176    if (!psVectorStats (stats, dT, NULL, NULL, 0)) {
     177        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     178        return false;
     179    }
    171180    float To = stats->robustMedian;
    172181    fprintf (stderr, "offset t: %f +/- %f\n", stats->robustMedian, stats->robustStdev);

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroLoadRefstars.c

    r23594 r24244  
    181181            ref->sky->d   = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC");
    182182            ref->Mag      = 0.001*psMetadataLookupS32 (&status, row, "MAG");  // ELIXIR uses millimags
     183            ref->Color    = 0.0;
    183184        } else {
    184185            ref->sky->r   = RAD_DEG*psMetadataLookupF64 (&status, row, "RA");
    185186            ref->sky->d   = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC");
    186187            ref->Mag      = psMetadataLookupF32 (&status, row, "MAG"); // PANSTARRS uses mags
     188            ref->Color    = 0.0;
    187189        }
    188190
    … …  
    219221        ref->sky->d   = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC");
    220222        ref->Mag      = psMetadataLookupF32 (&status, row, "MAG");
     223        float MagC1   = psMetadataLookupF32 (&status, row, "MAG_C1");
     224        float MagC2   = psMetadataLookupF32 (&status, row, "MAG_C2");
     225        if (!isnan(MagC1) && !isnan(MagC2)) {
     226            ref->Color = MagC1 - MagC2;
     227        } else {
     228            // XXX save the color and the slope in the table header?
     229            ref->Color = 0.0;
     230        }
    221231
    222232        // XXX VERY temporary hack to avoid M31 bulge

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroMaskUpdates.c

    r21409 r24244  
    1919}
    2020
    21 pmCell *pmCellInChip (pmChip *chip, float x, float y);
    22 bool pmCellCoordsForChip (float *xCell, float *yCell, pmCell *cell, float xChip, float yChip);
    23 bool pmChipCoordsForCell (float *xChip, float *yChip, pmCell *cell, float xCell, float yCell);
    24 
    25 bool psastroMaskCircle (psImage *mask, psImageMaskType value, float x0, float y0, float dX, float dY);
    26 bool psastroMaskBox (psImage *mask, psImageMaskType value, float x0, float y0, float dL, float dW, float theta);
    27 void psastroMaskLine (psImage *mask, psImageMaskType value, double x1, double y1, double x2, double y2, int dW);
    28 void psastroMaskLineBresen (psImage *mask, psImageMaskType value, int X1, int Y1, int X2, int Y2, int dW, int swapcoords);
    29 void psastroMaskRectangle (psImage *mask, psImageMaskType value, int x0, int y0, int x1, int y1);
    30 
    3121/**
    3222 * create a mask or mask regions based on the collection of reference stars that * are in the vicinity of each chip
    … …  
    4030    float zeropt, exptime;
    4131
    42     psImageMaskType maskValue  = pmConfigMaskGet("GHOST", config); // Mask value for ghost pixels
     32    psImageMaskType ghostMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for ghost pixels
     33    psImageMaskType spikeMaskValue = pmConfigMaskGet("SPIKE", config); // Mask value for ghost pixels
     34    psImageMaskType starMaskValue  = pmConfigMaskGet("STARCORE", config); // Mask value for ghost pixels
    4335
    4436    // psImageMaskType maskBlank  = pmConfigMaskGet("BLANK", config); // Mask value for blank pixels
    … …  
    5547    if (!REFSTAR_MASK) return true;
    5648
     49    // convert star positions to ghost positions and add to the readout->analysis data
     50    if (!psastroLoadGhosts (config)) {
     51        psError(PSASTRO_ERR_CONFIG, false, "Error loading ghosts");
     52        return false;
     53    }
     54    bool COUNT_GHOSTS = psMetadataLookupF32 (&status, recipe, "REFSTAR_COUNT_GHOSTS");
     55    double GHOST_MAX_MAG = psMetadataLookupF32 (&status, recipe, "GHOST_MAX_MAG");
     56    int nGhosts = 0;
     57
    5758    psLogMsg ("psastro", PS_LOG_INFO, "generating a bright-star mask");
    5859
    59     bool REFSTAR_MASK_REGIONS              = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_REGIONS");
     60    bool REFSTAR_MASK_BLEED                = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED");
    6061
    6162    double REFSTAR_MASK_MAX_MAG            = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_MAX_MAG");
    … …  
    9091    REFSTAR_MASK_SATSPIKE_MAG_MAX += MagOffset;
    9192    REFSTAR_MASK_BLEED_MAG_MAX += MagOffset;
     93    GHOST_MAX_MAG += MagOffset;
    9294
    9395    // select the output mask image :: we mosaic to chip mosaic format
    … …  
    132134
    133135        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE;
    134 
    135         // text region files for testing
    136         FILE *f = NULL;
    137         if (REFSTAR_MASK_REGIONS) {
    138           char *filename = NULL;
    139           char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME");
    140           psStringAppend (&filename, "refstars.mask.%s.dat", chipname);
    141           FILE *f = fopen (filename, "w");
    142           if (!f) {
    143             psWarning ("cannot create refstar mask file %s\n", filename);
    144             continue;
    145           }
    146           psFree (filename);
    147         }
    148136
    149137        pmChip *refChip  = pmFPAviewThisChip (view, refMask->fpa);
    … …  
    172160
    173161            // we mask pixels on the input mask image (chip-mosaic)
    174             pmCell *cellMask = pmFPAviewThisCell(view, outMask->fpa);
    175             pmReadout *readoutMask = NULL;
    176             if (cellMask->readouts->n) {
    177                 readoutMask = cellMask->readouts->data[0];
    178             }
     162            // pmCell *cellMask = pmFPAviewThisCell(view, outMask->fpa);
     163            // pmReadout *readoutMask = NULL;
     164            // if (cellMask->readouts->n) {
     165            //     readoutMask = cellMask->readouts->data[0];
     166            // }
    179167
    180168            // process each of the readouts
    … …  
    183171                if (! readout->data_exists) { continue; }
    184172
     173                // XXX why not do this?
     174                pmReadout *readoutMask = pmFPAviewThisReadout (view, outMask->fpa);
     175                if (!readoutMask) continue;
     176
    185177                // select the raw objects for this readout
    186178                psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
    187                 if (refstars == NULL) { continue; }
     179                if (!refstars) continue;
    188180
    189181                // we need to generate the following masks regions:
    … …  
    191183                // 2) diffraction spikes in direction ROT - ROTo
    192184                // 3) bleed trail in the direction of the readout
    193 
    194                 // XXX for the moment, let's just generate mana region-file objects in chip pixel space
     185                //    update: with 'burntool' applied to the data, the bleed trail mask is not needed.
     186
    195187                for (int i = 0; i < refstars->n; i++) {
    196188                    pmAstromObj *ref = refstars->data[i];
     189                   
     190                    if (COUNT_GHOSTS) {
     191                        if (ref->Mag > GHOST_MAX_MAG) {
     192                            nGhosts ++;
     193                        }
     194                    }
     195
    197196                    if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue;
    198197
    199                     // XXX convert ref->Mag to instrumental mags
     198                    // the reference magnitudes have been converted from the instrumental
     199                    // values supplied in the recipe to apparent mags (above)
    200200
    201201                    // CIRCLE around the stars (scaled by magnitude)
    202202                    float radius = REFSTAR_MASK_SATSTAR_MAG_SLOPE * (REFSTAR_MASK_SATSTAR_MAG_MAX - ref->Mag);
    203203
    204                     if (REFSTAR_MASK_REGIONS) {
    205                       fprintf (f, "CIRCLE %f %f  %f %f\n", ref->chip->x, ref->chip->y, radius, radius);
     204                    // XXX for now, assume cell binning is 1x1 relative to chip
     205                    psastroMaskCircle (readoutMask->mask, starMaskValue, ref->chip->x, ref->chip->y, radius, radius);
     206
     207                    for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) {
     208
     209                        float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO;
     210
     211                        // LINE for boundaries of the saturation spikes (scaled by magnitude)
     212                        // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2;
     213
     214                        float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX;
     215                        float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (MAG_MAX - ref->Mag);
     216                        float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH;
     217                        // XXX we can make the width depend on the spike as well...
     218                        // The length should also be a function of the image background level
     219
     220                        psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta);
    206221                    }
    207222
    208                     // XXX for now, assume cell binning is 1x1 relative to chip
    209                     if (readoutMask) {
    210                         psastroMaskCircle (readoutMask->mask, maskValue, ref->chip->x, ref->chip->y, radius, radius);
    211                     }
    212 
    213                     // LINE for boundaries of the saturation spikes (scaled by magnitude)
    214                     float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (REFSTAR_MASK_SATSPIKE_MAG_MAX - ref->Mag);
    215                     float spikeWidth = REFSTAR_MASK_SATSPIKE_WIDTH;
    216 
    217                     for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) {
    218                         float x0, y0, x1, y1, dx, dy;
    219 
    220                         float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO;
    221 
    222                         if (REFSTAR_MASK_REGIONS) {
    223                           // lower side
    224                           x0 = ref->chip->x + spikeWidth*sin(Theta);
    225                           y0 = ref->chip->y - spikeWidth*cos(Theta);
    226                           x1 = ref->chip->x + spikeLength*cos(Theta) + spikeWidth*sin(Theta);
    227                           y1 = ref->chip->y + spikeLength*sin(Theta) - spikeWidth*cos(Theta);
    228                           dx = x1 - x0;
    229                           dy = y1 - y0;
    230 
    231                           fprintf (f, "LINE %f %f  %f %f\n", x0, y0, dx, dy);
    232 
    233                           // upper side
    234                           x0 = ref->chip->x - spikeWidth*sin(Theta);
    235                           y0 = ref->chip->y + spikeWidth*cos(Theta);
    236                           x1 = ref->chip->x + spikeLength*cos(Theta) - spikeWidth*sin(Theta);
    237                           y1 = ref->chip->y + spikeLength*sin(Theta) + spikeWidth*cos(Theta);
    238                           dx = x1 - x0;
    239                           dy = y1 - y0;
    240                           fprintf (f, "LINE %f %f  %f %f\n", x0, y0, dx, dy);
    241                         }
    242 
    243                         if (readoutMask) {
    244                             psastroMaskBox (readoutMask->mask, maskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta);
    245                         }
    246                     }
    247 
    248                     // convert x,y chip coordinates to cells in maskChip
    249                     pmCell *refCell = pmCellInChip (refChip, ref->chip->x, ref->chip->y);
    250 
    251                     // LINE for boundaries of the bleed lines
    252                     if (REFSTAR_MASK_REGIONS) {
    253                       fprintf (f, "LINE %f %f  %f %f\n", ref->chip->x, ref->chip->y, 0.0, -100.0);
    254                     }
    255 
    256                     if (readoutMask && refCell) {
    257                         float xCell = 0.0;
    258                         float yCell = 0.0;
    259                         float xEnd = 0.0;
    260                         float yEnd = 0.0;
    261                         // find coordinate of star on cell
    262                         pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y);
    263                         // find coordinate of end-point on chip
    264 
    265                         int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE");
    266                         pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize);
    267 
    268                         float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
    269                         psastroMaskRectangle (readoutMask->mask, maskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd);
    270                     }
     223                    if (REFSTAR_MASK_BLEED) {
     224
     225                        // convert x,y chip coordinates to cells in maskChip
     226                        pmCell *refCell = pmCellInChip (refChip, ref->chip->x, ref->chip->y);
     227
     228                        // LINE for boundaries of the bleed lines
     229                        if (refCell) {
     230                            float xCell = 0.0;
     231                            float yCell = 0.0;
     232                            float xEnd = 0.0;
     233                            float yEnd = 0.0;
     234                            // find coordinate of star on cell
     235                            pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y);
     236                            // find coordinate of end-point on chip
     237
     238                            int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE");
     239                            pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize);
     240
     241                            float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
     242                            psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd);
     243                        }
     244                    }
    271245                }
     246
     247                // select the raw objects for this readout (loaded in psastroExtract.c)
     248                psArray *ghosts = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GHOSTS");
     249                if (ghosts == NULL) { continue; }
     250
     251                // identify the bright stars of interest
     252                for (int i = 0; i < ghosts->n; i++) {
     253                    psastroGhost *ghost = ghosts->data[i];
     254                    // XXX bright vs faint ghost bits? (OR with SUSPECT)
     255                    psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer);
     256                }
     257
     258                // select the raw objects for this readout, flag is they fall in a mask
     259                psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS");
     260                if (rawstars == NULL) return false;
     261               
     262                // XXX finish this:
     263                for (int i = 0; i < rawstars->n; i++) {
     264                    pmAstromObj *raw = rawstars->data[i];
     265                    psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[(int)(raw->chip->x)][(int)(raw->chip->y)];
     266                    if (value) continue;
     267                }
    272268            }
    273269        }
    274270
    275         // output sequence for mask corresponding to this chip
     271        // output sequence for mask corresponding to this chip (XXX this may not be needed...)
    276272        *viewMask = *view;
    277273        while ((cell = pmFPAviewNextCell (viewMask, outMask->fpa, 1)) != NULL) {
    … …  
    285281            if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE;
    286282        }
    287 
    288283        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
    289         if (REFSTAR_MASK_REGIONS) {
    290           fclose (f);
    291         }
    292284    }
    293285    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
     286
     287    if (COUNT_GHOSTS) {
     288        // save nGhosts to update header.
     289        psMetadata *updates = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
     290        if (!updates) {
     291            updates = psMetadataAlloc ();
     292            psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
     293            psFree (updates);
     294        }
     295        psMetadataAddS32 (updates, PS_LIST_TAIL, "NGHOSTS", PS_META_REPLACE, "total expected ghosts", nGhosts);
     296    }
    294297
    295298    // deactivate all files
    … …  
    301304}
    302305
    303 // XXX this is going to be very slow...
    304 pmCell *pmCellInChip (pmChip *chip, float x, float y) {
    305 
    306 # if (0)
    307     int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    308     int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    309     int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    310     int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
    311 
    312     // XXX fix the binning : currently not selected from concepts
    313     // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
    314     // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
    315     int xBin = 1;
    316     int yBin = 1;
    317 
    318     // Position on the cell
    319     float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin);
    320     float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin);
    321 # endif
    322 
    323     for (int i = 0; i < chip->cells->n; i++) {
    324 
    325         pmCell *cell = chip->cells->data[i];
    326         psRegion *region = pmCellExtent (cell);
    327 
    328         if (x < region->x0) goto skip;
    329         if (x > region->x1) goto skip;
    330         if (y < region->y0) goto skip;
    331         if (y > region->y1) goto skip;
    332 
    333         psFree (region);
    334         return cell;
    335 
    336     skip:
    337         psFree (region);
    338     }
    339     return NULL;
    340 }
    341 
    342 /**
    343  * convert chip coords to cell coords, given known cell
    344  */
    345 bool pmCellCoordsForChip (float *xCell, float *yCell, pmCell *cell, float xChip, float yChip) {
    346 
    347     int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    348     int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    349     int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    350     int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
    351 
    352     // XXX fix the binning : currently not selected from concepts
    353     // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
    354     // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
    355     int xBin = 1;
    356     int yBin = 1;
    357 
    358     // Position on the cell
    359     // ((pos)*(binning)*(cellParity) + (cell0))
    360     // XXX this is probably totally wrong now....
    361     // ((pos) - (cell0))*(cellParity)/(binning))
    362     *xCell = (xChip - x0Cell)*xParityCell/xBin;
    363     *yCell = (yChip - y0Cell)*yParityCell/yBin;
    364 
    365     return true;
    366 }
    367 
    368 /**
    369  * convert chip coords to cell coords, given known cell
    370  */
    371 bool pmChipCoordsForCell (float *xChip, float *yChip, pmCell *cell, float xCell, float yCell) {
    372 
    373     int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    374     int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    375     int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    376     int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
    377 
    378     // XXX fix the binning : currently not selected from concepts
    379     // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
    380     // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
    381     int xBin = 1;
    382     int yBin = 1;
    383 
    384     // Position on the cell
    385     // ((pos)*(binning)*(cellParity) + (cell0))
    386     // XXX this is probably totally wrong now....
    387     // ((pos) - (cell0))*(cellParity)/(binning))
    388     *xChip = xCell*xBin*xParityCell + x0Cell;
    389     *yChip = yCell*yBin*yParityCell + y0Cell;
    390 
    391     return true;
    392 }
    393 
    394 // XXX should be doing an OR
    395 bool psastroMaskCircle (psImage *mask, psImageMaskType value, float x0, float y0, float dX, float dY) {
    396 
    397     // XXX need to worry about row0, col0
    398     for (int ix = -dX; ix <= +dX; ix++) {
    399         int jx = ix + x0;
    400         if (jx < 0) continue;
    401         if (jx >= mask->numCols) continue;
    402         for (int iy = -dY; iy <= +dY; iy++) {
    403             int jy = iy + y0;
    404             if (jy < 0) continue;
    405             if (jy >= mask->numRows) continue;
    406 
    407             double r2 = PS_SQR(ix/dX) + PS_SQR(iy/dY);
    408             if (r2 > 1.0) continue;
    409 
    410             mask->data.PS_TYPE_IMAGE_MASK_DATA[jy][jx] |= value;
    411         }
    412     }
    413     return true;
    414 }
    415 
    416 // XXX should be doing an OR
    417 bool psastroMaskBox (psImage *mask, psImageMaskType value, float x0, float y0, float dL, float dW, float theta) {
    418 
    419     // draw a series of lines (from -0.5*dW to +0.5*dW) of length dL, starting at x0, y0, angle theta
    420 
    421     float xs = x0;
    422     float ys = y0;
    423 
    424     float xe = xs + dL*cos(theta);
    425     float ye = ys + dL*sin(theta);
    426 
    427     psastroMaskLine (mask, value, xs, ys, xe, ye, (int) dW);
    428     return true;
    429 }
    430 
    431 /**
    432  * identify the quadrant and draw the correct line
    433  */
    434 void psastroMaskLine (psImage *mask, psImageMaskType value, double x1, double y1, double x2, double y2, int dW) {
    435 
    436   int FlipDirect, FlipCoords;
    437   int X1, Y1, X2, Y2, dX, dY;
    438 
    439   /* rather than draw the line from float positions, we find the closest
    440      integer end-points and draw the line between those pixels */
    441 
    442   X1 = ROUND(x1);
    443   Y1 = ROUND(y1);
    444   X2 = ROUND(x2);
    445   Y2 = ROUND(y2);
    446 
    447   dX = X2 - X1;
    448   dY = Y2 - Y1;
    449 
    450   FlipCoords = (abs(dX) < abs(dY));
    451   FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2);
    452 
    453   if (!FlipDirect && !FlipCoords) psastroMaskLineBresen (mask, value, X1, Y1, X2, Y2, dW, FALSE);
    454   if ( FlipDirect && !FlipCoords) psastroMaskLineBresen (mask, value, X2, Y2, X1, Y1, dW, FALSE);
    455   if (!FlipDirect &&  FlipCoords) psastroMaskLineBresen (mask, value, Y1, X1, Y2, X2, dW, TRUE);
    456   if ( FlipDirect &&  FlipCoords) psastroMaskLineBresen (mask, value, Y2, X2, Y1, X1, dW, TRUE);
    457 
    458   return;
    459 }
    460 
    461 /**
    462  * use the Bresenham line drawing technique
    463  * integer-only Bresenham line-draw version which is fast
    464  */
    465 void psastroMaskLineBresen (psImage *mask, psImageMaskType value, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) {
    466 
    467     int X, Y, dX, dY;
    468     int e, e2;
    469 
    470     dX = X2 - X1;
    471     dY = Y2 - Y1;
    472 
    473     Y = Y1;
    474     e = 0;
    475     for (X = X1; X <= X2; X++) {
    476         if (X > 0) {
    477             if (swapcoords) {
    478                 if (X >= mask->numRows) continue;
    479                 for (int y = Y - dW; y <= Y + dW; y++) {
    480                     if (y < 0) continue;
    481                     if (y >= mask->numCols) continue;
    482                     mask->data.PS_TYPE_IMAGE_MASK_DATA[X][y] |= value;
    483                 }
    484             } else {
    485                 if (X >= mask->numCols) continue;
    486                 for (int y = Y - dW; y <= Y + dW; y++) {
    487                     if (y < 0) continue;
    488                     if (y >= mask->numRows) continue;
    489                     mask->data.PS_TYPE_IMAGE_MASK_DATA[y][X] |= value;
    490                 }
    491             }
    492         }
    493         e += dY;
    494         e2 = 2 * e;
    495         if (e2 > dX) {
    496             Y++;
    497             e -= dX;
    498         }
    499         if (e2 < -dX) {
    500             Y--;
    501             e += dX;
    502         }
    503     }
    504     return;
    505 }
    506 
    507 void psastroMaskRectangle (psImage *mask, psImageMaskType value, int x0, int y0, int x1, int y1) {
    508 
    509     int xs = PS_MAX (0, PS_MIN (mask->numCols, PS_MIN (x0, x1)));
    510     int xe = PS_MAX (0, PS_MIN (mask->numCols, PS_MAX (x0, x1)));
    511     int ys = PS_MAX (0, PS_MIN (mask->numRows, PS_MIN (y0, y1)));
    512     int ye = PS_MAX (0, PS_MIN (mask->numRows, PS_MAX (y0, y1)));
    513 
    514     for (int iy = ys; iy < ye; iy++) {
    515         for (int ix = xs; ix < xe; ix++) {
    516             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= value;
    517         }
    518     }
    519 }
    520 

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroMetadataStats.c

    r21409 r24244  
    1313# include "psastroInternal.h"
    1414
    15 bool psastroMetadataStats (pmConfig *config) {
     15bool psastroMetadataStats (pmConfig *config, psMetadata *stats) {
    1616
    1717    bool status;
    … …  
    2424
    2525    if (!output) {
    26         psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find output file (PSASTRO.OUTPUT).");
    27         return false;
     26        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find output file (PSASTRO.OUTPUT).");
     27        return false;
    2828    }
    2929
    30     // create output stats metadata
    31     psMetadata *stats = psMetadataAlloc ();
    32 
    33     // extract stats for the complete fpa
     30    // extract stats for the complete fpa
    3431    pmFPAview *view = pmFPAviewAlloc(0);
    3532
    3633    if (!ppStatsMetadata(stats, output->fpa, view, 0, config)) {
    37         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate stats for image.");
    38         psFree(view);
    39         psFree(stats);
    40         return false;
     34        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate stats for image.");
     35        psFree(view);
     36        return false;
    4137    }
    4238
    4339    // if we did not request any specific stats, the structure is empty
    4440    if (stats && stats->list->n == 0) {
    45         psWarning ("stats output specified, but no requested stats entries in headers");
    46         psFree(view);
    47         psFree(stats);
    48         return true;
     41        psWarning ("stats output specified, but no requested stats entries in headers");
     42        psFree(view);
     43        return true;
    4944    }
    5045
    … …  
    5247    char *statsMDC = psMetadataConfigFormat(stats);
    5348    if (!statsMDC || strlen(statsMDC) == 0) {
    54         psError(PS_ERR_IO, false, "Unable to serialize stats metadata.\n");
    55         return false;
    56     } 
     49        psError(PS_ERR_IO, false, "Unable to serialize stats metadata.\n");
     50        return false;
     51    }
    5752
    5853    // convert to a real UNIX filename
    … …  
    6055    FILE *statsFile = fopen (resolved, "w");
    6156    if (!statsFile) {
    62         psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
    63         psFree(stats);
    64         psFree(view);
    65         psFree(statsMDC);
    66         psFree(resolved);
    67         return false;
     57        psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
     58        psFree(view);
     59        psFree(statsMDC);
     60        psFree(resolved);
     61        return false;
    6862    }
    6963
    7064    // write the stats MDC to a file
    71     // XXX why does this not call psMetadataConfigPrint?
    7265    fprintf(statsFile, "%s", statsMDC);
    7366    fclose(statsFile);
     67
     68    pmConfigRunFilenameAddWrite(config, "STATS", filename);
    7469
    7570    psFree(resolved);
    7671    psFree(statsMDC);
    7772    psFree(view);
    78     psFree(stats);
    7973    return true;
    8074}

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroModelAnalysis.c

    r21409 r24244  
    153153
    154154    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    155     psVectorStats (stats, posZero, NULL, NULL, 0);
     155    if (!psVectorStats (stats, posZero, NULL, NULL, 0)) {
     156        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     157        return false;
     158    }
    156159
    157160    fprintf (outfile, "# pos zero %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroModelFitBoresite.c

    r21409 r24244  
    6262
    6363    // center (Xo) = mean(Xo), RX = range / 2
    64     psVectorStats (stats, Xo, NULL, NULL, 0);
     64    if (!psVectorStats (stats, Xo, NULL, NULL, 0)) {
     65        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     66        return NULL;
     67    }
    6568    params->data.F32[PAR_X0] = stats->sampleMean;
    6669    params->data.F32[PAR_RX] = (stats->max - stats->min) / 2.0;
    6770
    6871    // center (Yo) = mean(Yo), RY = range / 2
    69     psVectorStats (stats, Yo, NULL, NULL, 0);
     72    if (!psVectorStats (stats, Yo, NULL, NULL, 0)) {
     73        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     74        return NULL;
     75    }
    7076    params->data.F32[PAR_Y0] = stats->sampleMean;
    7177    params->data.F32[PAR_RY] = (stats->max - stats->min) / 2.0;

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroMosaicOneChip.c

    r21422 r24244  
    5858
    5959    // modify the order to correspond to the actual number of matched stars:
    60     if ((match->n < 17) && (order >= 3)) order = 2;
    61     if ((match->n < 13) && (order >= 2)) order = 1;
    62     if ((match->n <  9) && (order >= 1)) order = 0;
    63     if ((match->n <  3) || (order < 0) || (order > 3)) {
     60    int Ndof_min = 3;
     61    int order_max = 0.5*(sqrt(4*match->n - 4*Ndof_min + 1) - 3);
     62    order = PS_MIN (order, order_max);
     63
     64    // if ((match->n < 17) && (order >= 3)) order = 2;
     65    // if ((match->n < 13) && (order >= 2)) order = 1;
     66    // if ((match->n <  9) && (order >= 1)) order = 0;
     67
     68    if (order < 0) {
    6469        psLogMsg ("psastro", 3, "insufficient stars (%ld) or invalid order (%d)", match->n, order);
    6570        return false;

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroOneChipFit.c

    r21422 r24244  
    6464        // modify the order to correspond to the actual number of matched stars:
    6565        int Ndof_min = 3;
    66         int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1));
     66        int order_max = 0.5*(sqrt(4*match->n - 4*Ndof_min + 1) - 3);
    6767        order = PS_MIN (order, order_max);
     68
     69        // order 0 : Ro -> nterms = 1 * 2;
     70        // order 1 : Ro, Rx, Ry -> nterms = 3 * 2;
     71        // order 2 : Ro, Rx, Ry, Rxx, Rxy, Ryy -> nterms = 6 * 2;
     72        // order 3 : Ro, Rx, Ry, Rxx, Rxy, Ryy, Rxxx, Rxxy, Rxyy, Ryyy -> nterms = 10 * 2
     73        // 2*(N+1)*(N+2)/2 = (N+1)*(N+2) = nterms;
     74        // (order+1)(order+2) + ndof = nvalues
     75        // order^2 + 3*order + 2 + ndof = nvalue;
     76        // order^2 + 3*order + 2 + ndof - nvalue = 0;
     77        // 2*order = -3 +/- sqrt (9 - 4*(2 - nvalue + ndof));
     78        // 2*order = -3 +/- sqrt (9 - 8 + 4*nvalue - 4*ndof);
     79        // 2*order = (sqrt (1 + 4*nvalue - 4*ndof) - 3);
    6880
    6981        // if ((match->n < 11) && (order >= 3)) order = 2;

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroParseCamera.c

    r21409 r24244  
    4545    psFree (chips);
    4646
    47     psString dump_file = psMetadataLookupStr(&status, config->arguments, "DUMP_CONFIG");
    48     if (dump_file) {
    49         pmConfigCamerasCull(config, NULL);
    50         pmConfigRecipesCull(config, "PPIMAGE,PPSTATS,PSPHOT,MASKS,PSASTRO");
    51 
    52         pmConfigDump(config, input->fpa, dump_file);
    53     }
    54 
    55 
    5647    psTrace("psastro", 1, "Done with psastroParseCamera...\n");
    5748    return true;

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroVersion.c

    r23352 r24244  
    99
    1010#include "psastroInternal.h"
     11#include "psastroVersionDefinitions.h"
    1112
    1213#ifndef PSASTRO_VERSION
    … …  
    2021#endif
    2122
    22 #define xstr(s) str(s)
    23 #define str(s) #s
    24 
    2523psString psastroVersion(void)
    2624{
    2725    char *value = NULL;
    28     psStringAppend(&value, "%s@%s", xstr(PSASTRO_BRANCH), xstr(PSASTRO_VERSION));
     26    psStringAppend(&value, "%s@%s", PSASTRO_BRANCH, PSASTRO_VERSION);
    2927    return value;
    3028}
    … …  
    3230psString psastroSource(void)
    3331{
    34   return psStringCopy(xstr(PSASTRO_SOURCE));
     32  return psStringCopy(PSASTRO_SOURCE);
    3533}
    3634

  • branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroZeroPoint.c

    r21409 r24244  
    126126    // this analysis has too few data points to use the robust median method
    127127    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    128     psVectorStats (stats, dMag, NULL, NULL, 0);
     128    if (!psVectorStats (stats, dMag, NULL, NULL, 0)) {
     129        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     130        return false;
     131    }
    129132    fprintf (stderr, "zero point %f +/- %f using %d stars; transparency %f\n", stats->clippedMean, stats->clippedStdev, Npts, zeropt - stats->clippedMean);
    130133
Note: See TracChangeset for help on using the changeset viewer.