IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15363


Ignore:
Timestamp:
Oct 23, 2007, 11:50:52 AM (19 years ago)
Author:
bills
Message:

Extract postage stamp from chip mosaic image. Changed -celcenter to -skycenter
Range in celestial coordinates not specified with -arcrange (arcseconds)
Better error handling.
Allow files without WCS if we don't need it

Location:
trunk/ppstamp/src
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppstamp/src

    • Property svn:ignore

      •  

        old new  
        66config.h.in
        77stamp-h1
         8ppstampErrorCodes.c
         9ppstampErrorCodes.h

  • trunk/ppstamp/src/.cvsignore

    r15280 r15363  
    66config.h.in
    77stamp-h1
     8ppstampErrorCodes.c
     9ppstampErrorCodes.h

  • trunk/ppstamp/src/Makefile.am

    r15323 r15363  
    1818        ppstampErrorCodes.c \
    1919        ppstampMakeStamp.c \
     20        ppstampMosaic.c \
    2021        ppstampOptions.c \
    2122        ppstampParseCamera.c \

  • trunk/ppstamp/src/ppstamp.c

    r15280 r15363  
    77int main(int argc, char **argv)
    88{
    9     ppstampOptions *options;
     9    ppstampOptions *options = NULL;
     10    int exitCode;
    1011
    1112    psTimerStart(TIMERNAME);
    12 
    1313    psLibInit(NULL);
    14 
    15     // Parse the configuration and arguments
    16     // Open the input image(s)
    17     // note usage is called on error or request by ppstampArguments
    18     // TODO is the cleanup really necessary since we're exiting anyways?
    19     // (It's useful to know about memory leaks)
     14    ppstampErrorRegister();
    2015
    2116    pmConfig *config = ppstampArguments(argc, argv, &options);
    … …  
    2621    }
    2722
    28     // define recipe options
    2923    // define the active I/O files
    3024    if (!ppstampParseCamera(config)) {
    … …  
    3529
    3630    // find the pixels that we need to copy, setup the output image
    37     if (!ppstampMakeStamp(config, options)) {
    38         exit(PS_EXIT_UNKNOWN_ERROR);
     31    if (ppstampMakeStamp(config, options)) {
     32        exitCode = 0;
     33    } else {
     34        exitCode = PS_EXIT_DATA_ERROR;
    3935    }
    4036
    … …  
    4440    ppstampCleanup(config, options);
    4541
    46     return PS_EXIT_SUCCESS;
     42    return exitCode;
    4743}

  • trunk/ppstamp/src/ppstamp.h

    r15323 r15363  
    1515#include "psastro.h"
    1616#include "ppStats.h"
     17#include "ppstampErrorCodes.h"
    1718
    1819#define TIMERNAME "ppstamp"             // Name of timer
    … …  
    5960
    6061bool ppstampMakeStamp(pmConfig *config, ppstampOptions *);
     62pmFPAfile * ppstampBuildMosaic(pmConfig *config, pmFPAfile *input, pmFPAview *view);
    6163
    6264// free memory, check for leaks
    … …  
    7072
    7173/// Update the metadata with version information for all dependencies
    72 void ppstampVersionMetadata(psMetadata *metadata ///< Metadata to update with version information
    73     );
     74void ppstampVersionMetadata(psMetadata *metadata, ppstampOptions *options);
    7475
    7576psRegion *ppstampCellRegion(const pmCell *cell);

  • trunk/ppstamp/src/ppstampArguments.c

    r15323 r15363  
    66#include "ohana.h"
    77
    8 static bool get2Angles(int argnum, int *pArgc, char **argv, bool bothDegrees, double *p1, double *p2);
    9 static bool get2Ints(int argnum, int *pArgc, char **argv, int *p1, int *p2);
    10 
    11 // XXX This usage output is kind of busy
     8static bool get2Angles(int argnum, int *pArgc, char **argv, bool bothDegrees, bool makePositive, double *p1, double *p2);
     9static bool get2Ints(int argnum, int *pArgc, char **argv, bool makePositive, int *p1, int *p2);
    1210
    1311static void usage (void)
    … …  
    1513    fprintf(stderr, "\n");
    1614
    17     fprintf(stderr, "USAGE: ppstamp -celcenter RA DEC -arcrange dRA dDEC [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
    18     fprintf(stderr, "       ppstamp -celcenter RA DEC -celrange dRA dDEC [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
    19     fprintf(stderr, "       ppstamp -celcenter RA DEC -pixrange dx dy    [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
    20     fprintf(stderr, "       ppstamp -pixcenter x y    -pixrange dx dy    -chip chipname [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
    21     fprintf(stderr, "       ppstamp -pixcenter x y    -arcrange dRA dDEC -chip chipname [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
    22     fprintf(stderr, "       ppstamp -pixcenter x y    -celrange dRA dDEC -chip chipname [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
    23     fprintf(stderr, "\n");
    24 
    25     fprintf(stderr, "celcenter & celrange may be specified in HH:MM:SS DD:MM:SS or as 2 values in decimal degrees\n");
    26     fprintf(stderr, "arcrange may be specfied as 2 values DD:MM:SS DD:MM:SS or decimal degrees\n");
     15    fprintf(stderr, "USAGE: ppstamp -skycenter RA DEC -pixrange dx dy    [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
     16    fprintf(stderr, "       ppstamp -skycenter RA DEC -arcrange dRA dDEC [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
     17    fprintf(stderr, "       ppstamp -pixcenter x y    -pixrange dx dy    [-chip chipname] [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
     18    fprintf(stderr, "       ppstamp -pixcenter x y    -arcrange dRA dDEC [-chip chipname] [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
     19    fprintf(stderr, "\n");
     20
     21    fprintf(stderr, "-skycenter is specified as HH:MM:SS DD:MM:SS or decimal degrees\n");
     22    fprintf(stderr, "-arcrange values are seconds of arc\n");
    2723    fprintf(stderr, "\n");
    2824    fprintf(stderr, "Optional arguments:\n");
    … …  
    5551    pmConfig *config = pmConfigRead(&argc, argv, NULL);
    5652    if (config == NULL) {
    57         psErrorStackPrint(stderr, "Can't find site configuration!\n");
    58         exit(EXIT_FAILURE);
     53        psError(PPSTAMP_ERR_CONFIG, false, "Can't read site configuration!\n");
     54        return(NULL);
    5955    }
    6056
    … …  
    6763        psArgumentRemove(argnum, &argc, argv);
    6864
    69         if (!get2Ints(argnum, &argc, argv, &options->centerX, &options->centerY)) {
    70             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "invalid pixcenter specification\n");
    71             usage();
    72         }
    73     }
    74 
    75     if ((argnum = psArgumentGet(argc, argv, "-celcenter"))) {
     65        if (!get2Ints(argnum, &argc, argv, false, &options->centerX, &options->centerY)) {
     66            psError(PPSTAMP_ERR_ARGUMENTS, true, "invalid pixcenter specification: %s %s\n",
     67                argv[argnum], argv[argnum+1]);
     68            usage();
     69        }
     70        psArgumentRemove(argnum, &argc, argv);
     71        psArgumentRemove(argnum, &argc, argv);
     72    }
     73
     74    if ((argnum = psArgumentGet(argc, argv, "-skycenter"))) {
    7675        if (gotCenter) {
    77             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "can't specify both -pixcenter and -celcenter\n");;
     76            psError(PPSTAMP_ERR_ARGUMENTS, true, "can't specify both -pixcenter and -skycenter\n");;
    7877            usage();
    7978        }
    … …  
    8180
    8281        double raDeg, decDeg;
    83         if (!get2Angles(argnum, &argc, argv, false, &raDeg, &decDeg)) {
    84             usage();
    85         }
     82        if (!get2Angles(argnum, &argc, argv, false, false, &raDeg, &decDeg)) {
     83            psError(PPSTAMP_ERR_ARGUMENTS, true, "invalid skycenter specification: %s %s\n",
     84                argv[argnum], argv[argnum+1]);
     85            usage();
     86        }
     87        psArgumentRemove(argnum, &argc, argv);
     88        psArgumentRemove(argnum, &argc, argv);
    8689        options->centerRA  = DEG_TO_RAD(raDeg);
    8790        options->centerDEC = DEG_TO_RAD(decDeg);
    … …  
    9194
    9295    if (!gotCenter) {
    93         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify center of region of interest\n");
     96        psError(PPSTAMP_ERR_ARGUMENTS, true, "must specify center of region of interest\n");
    9497        usage();
    9598    }
    … …  
    99102        options->celestialRange = false;
    100103        psArgumentRemove(argnum, &argc, argv);
    101         if (!get2Ints(argnum, &argc, argv, &options->dX, &options->dY)) {
    102             usage();
    103         }
     104        if (!get2Ints(argnum, &argc, argv, true, &options->dX, &options->dY)) {
     105            psError(PPSTAMP_ERR_ARGUMENTS, true, "invalid pixrange specification: %s %s\n",
     106                argv[argnum], argv[argnum+1]);
     107            usage();
     108        }
     109        psArgumentRemove(argnum, &argc, argv);
     110        psArgumentRemove(argnum, &argc, argv);
    104111    }
    105112
    106113    if ((argnum = psArgumentGet(argc, argv, "-arcrange"))) {
    107114        if (gotRange) {
    108             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "specify only one of -pixrange or -arcrange or -celrange\n");;
     115            psError(PPSTAMP_ERR_ARGUMENTS, true, "specify only one of -pixrange or -arcrange\n");;
    109116            usage();
    110117        }
    … …  
    113120        options->celestialRange = true;
    114121
     122        // arcrange values are seconds of arc
     123        options->dRA  = DEG_TO_RAD(atof(argv[argnum]) * 3600.);
     124        options->dDEC = DEG_TO_RAD(atof(argv[argnum+1]) * 3600.);
     125
     126        psArgumentRemove(argnum, &argc, argv);
     127        psArgumentRemove(argnum, &argc, argv);
     128    }
     129    // I'm leaving in the -celrange option (HH:MM:SS DD:MM:SS), but not publicizing it
     130    if ((argnum = psArgumentGet(argc, argv, "-celrange"))) {
     131        if (gotRange) {
     132            psError(PPSTAMP_ERR_ARGUMENTS, true, "specify one of -pixrange, -arcrange, or -celrange\n");;
     133            usage();
     134        }
     135        psArgumentRemove(argnum, &argc, argv);
     136        gotRange = true;
     137        options->celestialRange = true;
     138
    115139        double deg1, deg2;
    116         if (!get2Angles(argnum, &argc, argv, true, &deg1, &deg2)) {
    117             usage();
    118         }
     140        if (!get2Angles(argnum, &argc, argv, false, true, &deg1, &deg2)) {
     141            psError(PPSTAMP_ERR_ARGUMENTS, true, "invalid celrange specification: %s %s\n",
     142                argv[argnum], argv[argnum+1]);
     143            usage();
     144        }
     145        psArgumentRemove(argnum, &argc, argv);
     146        psArgumentRemove(argnum, &argc, argv);
    119147        options->dRA = DEG_TO_RAD(deg1);
    120148        options->dDEC = DEG_TO_RAD(deg2);
    121149    }
    122     if ((argnum = psArgumentGet(argc, argv, "-celrange"))) {
    123         if (gotRange) {
    124             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "specify one of -pixrange or -arcrange or -celrange\n");;
    125             usage();
    126         }
    127         psArgumentRemove(argnum, &argc, argv);
    128         gotRange = true;
    129         options->celestialRange = true;
    130 
    131         double deg1, deg2;
    132         if (!get2Angles(argnum, &argc, argv, false, &deg1, &deg2)) {
    133             usage();
    134         }
    135         options->dRA = DEG_TO_RAD(deg1);
    136         options->dDEC = DEG_TO_RAD(deg2);
    137     }
    138150
    139151    if (!gotRange) {
    140         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify range\n");
     152        psError(PPSTAMP_ERR_ARGUMENTS, true, "must specify range\n");
    141153        usage();
    142154    }
    … …  
    151163    bool status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list");
    152164    if (!status) {
    153         // TODO: shall we print a specific message?
    154         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify INPUT\n");
    155         usage ();
     165        psError(PPSTAMP_ERR_ARGUMENTS, true, "must specify INPUT\n");
     166        usage();
    156167    }
    157168
    158169    // finally the output file
    159     if (argc != 2) {
    160         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify OUTPUT\n");
    161         usage ();
    162     }
    163 
    164     // Add the input and output images (which remain on the command-line) to the arguments list
     170    if (argc < 2) {
     171        psError(PPSTAMP_ERR_ARGUMENTS, true, "must specify OUTPUT\n");
     172        usage();
     173    } else if (argc > 2) {
     174        psError(PPSTAMP_ERR_ARGUMENTS, true, "unknown argument(s)\n");
     175        usage();
     176    }
     177
     178    // Add the output image (which remains on the command-line) to the arguments list
    165179    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]);
    166180
    … …  
    173187{
    174188    char *p = string;
     189
    175190    if ((*p == '+') || (*p == '-')) {
    176191        p++;
    … …  
    179194}
    180195
    181 static bool get2Ints(int argnum, int *pArgc, char **argv, int *p1, int *p2)
     196static bool get2Ints(int argnum, int *pArgc, char **argv, bool makePositive, int *p1, int *p2)
    182197{
    183198    if (*pArgc < 2) {
    … …  
    188203        return false;
    189204    }
    190 
    191205    *p1 = atoi(argv[argnum]);
    192     psArgumentRemove(argnum, pArgc, argv);
    193 
    194     if (!validNumber(argv[argnum])) {
     206
     207    if (!validNumber(argv[argnum+1])) {
    195208        return false;
    196209    }
    197 
    198     *p2 = atoi(argv[argnum]);
    199     psArgumentRemove(argnum, pArgc, argv);
     210    *p2 = atoi(argv[argnum+1]);
     211
     212    if (makePositive) {
     213        *p1 = abs(*p1);
     214        *p2 = abs(*p2);
     215    }
    200216
    201217    return true;
    202218}
    203219
    204 static bool get2Angles(int argnum, int *pArgc, char **argv, bool bothInDegrees, double *p1, double *p2)
     220static bool get2Angles(int argnum, int *pArgc, char **argv, bool bothInDegrees, bool makePositive,
     221    double *p1, double *p2)
    205222{
    206223    bool rval;
    … …  
    221238    }
    222239
    223     psArgumentRemove(argnum, pArgc, argv);
    224     psArgumentRemove(argnum, pArgc, argv);
     240    if (rval && makePositive) {
     241        *p1 = abs(*p1);
     242        *p2 = abs(*p2);
     243    }
    225244
    226245    return rval;

  • trunk/ppstamp/src/ppstampCleanup.c

    r15318 r15363  
    2020
    2121    // psMemCheckLeaks (0, NULL, stderr, false);
    22     fprintf (stderr, "Found %d leaks in %s\n", psMemCheckLeaks (0, NULL, NULL, false), "ppstamp");
     22    // fprintf (stderr, "Found %d leaks in %s\n", psMemCheckLeaks (0, NULL, NULL, false), "ppstamp");
    2323
    2424    return;

  • trunk/ppstamp/src/ppstampErrorCodes.h.in

    r15280 r15363  
    99 */
    1010typedef enum {
    11     PPSTAMP_ERR_BASE = 600,
     11    PPSTAMP_ERR_BASE = 650,
    1212    PPSTAMP_ERR_${ErrorCode},
    1313    PPSTAMP_ERR_NERROR
    14 } pswarpErrorCode;
     14} ppstampErrorCode;
    1515
    16 void pswarpErrorRegister(void);
     16void ppstampErrorRegister(void);
    1717
    1818#endif

  • trunk/ppstamp/src/ppstampMakeStamp.c

    r15323 r15363  
    1616} ppstampOverlap;
    1717
     18static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);
     19
    1820// convert the input chip's transforms to the output
    1921static bool convertWCS(pmHDU *outHDU, pmFPA *outFPA, pmChip *outChip, pmChip *inChip, psRegion *roi)
    … …  
    2426
    2527    outChip->toFPA   = psPlaneTransformSetCenter(NULL, inChip->toFPA, roi->x0, roi->y0);
     28
    2629    outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50);
    2730
    2831    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_NONLIN_TOL)) {
    29         psError(PS_ERR_UNKNOWN, false, "Failed to write WCS\n");
     32        psError(PS_ERR_UNKNOWN, false, "Failed to write WCS to output\n");
    3033        return false;
    3134    }
    … …  
    3437}
    3538
    36 static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, pmCell *inCell, psRegion *roi)
     39static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, ppstampOptions *options)
    3740{
    3841    pmChip    *outChip;
    … …  
    4447    psFree(view);
    4548
    46     // XXX we probably should copy some concepts from the chip (skipping those that don't apply)
     49    // copy data from the input chip header to the output.
     50    // since some of the keywords might be duplicated we may not want to copy both
     51
     52    // copy the fpa concepts
    4753    outChip->parent->concepts = psMetadataCopy(outChip->parent->concepts, inChip->parent->concepts);
    4854
    49     pmHDU *inHDU  = pmHDUGetHighest(input->fpa, inChip, NULL);
     55    pmHDU *inHDU  = pmHDUFromChip(inChip);
    5056    pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);
    5157
    … …  
    5662    }
    5763
    58     // steal the input fpa's transforms
    59     output->fpa->toTPA = input->fpa->toTPA;
    60     output->fpa->fromTPA = input->fpa->fromTPA;
    61     output->fpa->toSky = input->fpa->toSky;
    62 
    63     // drop the references
    64     input->fpa->toTPA = 0;
    65     input->fpa->fromTPA = 0;
    66     input->fpa->toSky = 0;
    67 
    68     if (!convertWCS(outHDU, output->fpa, outChip, inChip, roi)) {
    69         return false;
    70     }
     64    // If input had WCS convert it for stamp
     65    if (input->fpa->toSky) {
     66        // steal the input fpa's transforms
     67        output->fpa->toTPA   = input->fpa->toTPA;
     68        output->fpa->fromTPA = input->fpa->fromTPA;
     69        output->fpa->toSky   = input->fpa->toSky;
     70
     71        // drop the references
     72        input->fpa->toTPA   = 0;
     73        input->fpa->fromTPA = 0;
     74        input->fpa->toSky   = 0;
     75
     76        if (!convertWCS(outHDU, output->fpa, outChip, inChip, &options->roi)) {
     77            return false;
     78        }
     79    } else {
     80        psWarning("No WCS present in input\n");
     81    }
     82
     83    ppstampVersionMetadata(outHDU->header, options);
    7184
    7285    return true;
    7386}
    7487
     88#ifdef notdef
     89// update the ROI to account for any change in the coordinate system in the
     90// mosaic process
     91// THERE isn't any change. This is something I was trying to make the megacam
     92// transforms come out correctly.
     93static bool updateROI(ppstampOptions *options, pmFPAfile *mosaic, pmChip *mosaicChip)
     94{
     95    // set up the astrometry
     96    pmHDU *hdu = pmHDUFromChip(mosaicChip);
     97    PS_ASSERT_PTR_NON_NULL(hdu, 1)
     98    PS_ASSERT_PTR_NON_NULL(hdu->header, 1)
     99
     100
     101    if (!pmAstromReadWCS(mosaic->fpa, mosaicChip, hdu->header, 1.0)) {
     102        psError(PS_ERR_UNKNOWN, false, "Failed to Read WCS from mosaic, cannot proceed\n");
     103        return false;
     104    }
     105    pmAstromObj *center = pmAstromObjAlloc();
     106
     107    center->sky->r = options->centerRA;
     108    center->sky->d = options->centerDEC;
     109    center->sky->rErr = 0;
     110    center->sky->dErr = 0;
     111
     112    skyToChip(center, mosaic->fpa, mosaicChip);
     113
     114    int width = options->dX;
     115    int height = options->dY;
     116
     117    fprintf(stderr, "ROI before: %s\n", psRegionToString(options->roi));
     118
     119    options->roi.x0 = center->chip->x - width/2;
     120    options->roi.x1 = options->roi.x0 + width;
     121    options->roi.y0 = center->chip->y - height/2;
     122    options->roi.y1 = options->roi.y0 + height;
     123
     124    fprintf(stderr, "ROI after:  %s\n", psRegionToString(options->roi));
     125
     126
     127    return true;
     128}
     129#endif
     130
    75131// Build the postage stamp output file
    76132
    77133static bool makeStamp(pmConfig *config, ppstampOptions *options, pmFPAfile *input,
    78                 pmChip *inChip, pmCell *inCell, pmFPAview *view)
     134                pmChip *inChip, pmFPAview *view)
    79135{
    80136    int status = false;
    … …  
    91147    outview->chip = 0;
    92148    outview->cell = 0;
    93     outview->readout = 0;
    94149    pmCell *target =  pmFPAviewThisCell(outview, output->fpa); // Target cell
    95150    psFree(outview);
    96151    outview = NULL;
    97152
    98     //    psMetadataPrint(stderr, inChip->concepts, 0);
    99     //   psMetadataPrint(stderr, inCell->concepts, 0);
    100     // these default to zero would that be ok?
     153    //   psMetadataPrint(stderr, inChip->concepts, 0);
     154   
     155    // These default to zero. would that be ok?
    101156    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_META_REPLACE, "Binning in x", 1);
    102157    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_META_REPLACE, "Binning in y", 1);
    103158
     159
     160    // we extract our postage stamp from a mosaic so that
     161    // pmChipMosaic can handle the tricky bits of parity, binning, etc.
     162    pmFPAfile *mosaic = ppstampBuildMosaic(config, input, view);
     163    if (mosaic == NULL) {
     164        return false;
     165    }
     166
     167    // one cell one chip
     168    pmFPAview *mosaicView = pmFPAviewAlloc(0);
     169    mosaicView->chip = 0;
     170    mosaicView->cell = 0;
     171
    104172    pmReadout *readout;
    105     while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
     173    while ((readout = pmFPAviewNextReadout (mosaicView, mosaic->fpa, 1)) != NULL) {
    106174        if (!readout->data_exists) {
     175            psError(PS_ERR_UNKNOWN, false, "no data in mosaic readout!\n");
    107176            continue;
    108177        }
    … …  
    113182            break;
    114183        }
    115 
     184   
    116185        psImage *subsetImage = psImageSubset(readout->image, options->roi);
    117186        if (subsetImage) {
    … …  
    120189            psFree(subsetImage);
    121190            if (outReadout->image) {
    122                 // The image has inherited the source's col0 and row0. We don't want this in
    123                 // our output so override it.
     191                // The image has inherited the subset's col0 and row0. We don't want this in
     192                // our output so override the values.
    124193                // XXX put this into a function in psImage
    125                 // (What I really needed was a function to "create a copy of this portion of an image")
     194                // (What I really needed was a function that does
     195                //     "create a copy of this portion of an image and return it as a new image with
     196                //      origin 0,0 ")
    126197                P_PSIMAGE_SET_COL0(outReadout->image, 0);
    127198                P_PSIMAGE_SET_ROW0(outReadout->image, 0);
    128199
    129200                outReadout->data_exists = true;
    130                 outReadout->parent->data_exists = true;         // cell
    131                 outReadout->parent->parent->data_exists = true; // chip
     201                outReadout->parent->data_exists = true;
     202                outReadout->parent->parent->data_exists = true;
    132203                status = true;
    133204            } else {
    … …  
    142213    }
    143214
     215    // XXXX remove this or make writing the mosaic to a file a valid option
     216    static bool writeMosaic = false;
     217    if (writeMosaic) {
     218        mosaic->save = true;
     219        mosaicView->cell = -1;
     220        pmFPAfileIOChecks(config, mosaicView, PM_FPA_AFTER);
     221    }
     222
     223    psFree(mosaicView);
     224
    144225    if (status) {
    145         status = copyMetadata(output, input, inChip, inCell, &options->roi);
     226        status = copyMetadata(output, input, inChip, options);
    146227    }
    147228    return status;
    … …  
    181262static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
    182263{
    183     // convert from sky to FP
     264    // convert from sky to TP to FP
    184265    psProject(pt->TP, pt->sky, fpa->toSky);
    185266    psPlaneTransformApply(pt->FP, fpa->fromTPA, pt->TP);
    … …  
    190271static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
    191272{
     273    // chip to FP
    192274    psPlaneTransformApply(pt->FP, chip->toFPA, pt->chip);
     275    // FP to TP to sky
    193276    psPlaneTransformApply(pt->TP, fpa->toTPA, pt->FP);
    194277    psDeproject(pt->sky, pt->TP, fpa->toSky);
    … …  
    208291}
    209292
     293// For roi width and height given in sky coordinates find a bounding box in chip coordinates
     294// that encloses the requested range
    210295static void findBoundingBox(ppstampOptions *options, pmFPA *fpa, pmChip *chip, pmAstromObj *center)
    211296{
    … …  
    249334
    250335    psFree(pt);
    251 }
    252 
    253 
     336
     337    // save the width and height in case we need them later
     338    options->dX = options->roi.x1 - options->roi.x0;
     339    options->dY = options->roi.y1 - options->roi.y0;
     340}
     341
     342#ifdef notdef
     343// I'm not using this any more. We mosaic the chip before extracting the pixels
     344// so we don't need to deal with cells at this level.
     345
     346// findCell
    254347// find cell on given chip that contains the region of interest
    255348// return the status of the overlap and if the cell completely contains the ROI
    … …  
    263356    view->cell = -1;
    264357    while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
    265         psRegion *cellExtent = ppstampCellRegion(cell);
     358        psRegion *cellExtent = pmCellExtent(cell);
    266359        if (regionContainsPoint(cellExtent, center->chip)) {
    267360            if (regionContainsRegion(cellExtent, &options->roi)) {
    … …  
    288381    return PPSTAMP_OFF;
    289382}
     383#endif
     384
     385
    290386
    291387// findROI
    292388// calculate the region of interest in chip coordinates and determine whether that region
    293 // whether the region of interest is completely contained in a cell on a given chip.
     389// is completely contained on a single chip
    294390
    295391static ppstampOverlap findROI(ppstampOptions *options, pmFPAview *view, pmFPAfile *input, pmChip *chip,
    296                        pmAstromObj *center, pmCell **ppCell)
     392                       pmAstromObj *center)
    297393{
    298394    psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
    299     psRegion    *chipExtent = ppstampChipRegion(chip);
     395    psRegion    *chipBounds = ppstampChipRegion(chip);
    300396    bool        onChip = false;
    301397    ppstampOverlap   returnval = PPSTAMP_OFF;
     398
     399//    fprintf(stderr, "ppstampChipBounds: %s\n", psRegionToString(*chipBounds));
    302400
    303401    // set up the astrometry
    … …  
    307405
    308406    if (!pmAstromReadWCS(input->fpa, chip, hdu->header, 1.0)) {
    309         psError(PS_ERR_UNKNOWN, false, "Failed to Read WCS\n");
    310         psFree(chipExtent);
    311         return PPSTAMP_ERROR;
     407        // we can live without WCS if the ROI is specified in pixels
     408
     409        if (options->celestialCenter || options->celestialRange) {
     410            psError(PPSTAMP_ERR_DATA, false, "Failed to Read WCS, cannot proceed\n");
     411            psFree(chipBounds);
     412            return PPSTAMP_ERROR;
     413        }
    312414    }
    313415
    … …  
    320422
    321423        skyToChip(center, input->fpa, chip);
    322         psPlaneTransformApply(center->chip, chip->fromFPA, center->FP);
    323 
    324         if (regionContainsPoint(chipExtent, center->chip)) {
    325             psLogMsg("ppstampMakeStamp", 3, "Found center (%f %f) on chip: %s\n",
     424
     425        if (regionContainsPoint(chipBounds, center->chip)) {
     426            psLogMsg("ppstampMakeStamp", 2, "Found center (%f %f) on chip: %s\n",
    326427                center->chip->x, center->chip->y, chipName);
    327428            onChip = true;
    … …  
    332433        // If no chip name was specified, select this one (the first one that had data)
    333434        if ((options->chipName == NULL) || !strcmp(chipName, options->chipName)) {
    334             psLogMsg("ppstampMakeStamp", 3, "Center on chip: %s\n", chipName);
     435            psLogMsg("ppstampMakeStamp", 2, "Center on chip: %s\n", chipName);
    335436            center->chip->x = options->centerX;
    336437            center->chip->y = options->centerY;
    … …  
    355456        }
    356457
    357         if (regionContainsRegion(chipExtent, &options->roi)) {
    358             returnval = findCell(view, options, input, center, ppCell);
     458
     459        if (regionContainsRegion(chipBounds, &options->roi)) {
     460            // returnval = findCell(view, options, input, center, ppCell);
     461            psLogMsg("ppstampMakeStamp", 2, "ROI contained on: %s\n", chipName);
     462            returnval = PPSTAMP_ON;
    359463        } else {
    360464            fprintf(stderr, "Partial Overlap chip %s\n", chipName);
    361465            fprintf(stderr, "ROI:         %s\n", psRegionToString(options->roi));
    362             fprintf(stderr, "Chip Extent: %s\n", psRegionToString(*chipExtent));
     466            fprintf(stderr, "Chip Extent: %s\n", psRegionToString(*chipBounds));
    363467
    364468            returnval = PPSTAMP_PARTIALLY_ON;
    365         }
    366     }
    367     psFree(chipExtent);
     469
     470        }
     471    }
     472    psFree(chipBounds);
    368473
    369474    return returnval;
    … …  
    374479    bool        status = false;
    375480    bool        returnval = false;;
     481    bool        foundOverlap = false;
    376482    pmAstromObj *center = pmAstromObjAlloc();
    377483
    … …  
    408514        }
    409515
    410         pmCell *cell;
    411         ppstampOverlap overlap = findROI(options, view, input, chip, center, &cell);
    412 
    413         if (overlap == PPSTAMP_ON) {
    414 
    415             returnval = makeStamp(config, options, input, chip, cell, view);
    416 
     516        ppstampOverlap overlap = findROI(options, view, input, chip, center);
     517
     518        switch (overlap) {
     519        case PPSTAMP_OFF:
     520            // keep looking
     521            break;
     522        case PPSTAMP_ON:
     523            returnval = makeStamp(config, options, input, chip, view);
    417524            allDone = true;
    418         } else if (overlap == PPSTAMP_PARTIALLY_ON) {
    419             psError(PS_ERR_UNKNOWN, false, "Region of interest must be confined to a single cell\n");
     525            foundOverlap = true;
     526            break;
     527        case PPSTAMP_PARTIALLY_ON:
     528            psError(PS_ERR_UNKNOWN, false, "Region of interest must be confined to a single chip\n");
    420529            // XXX: perhaps print a helpful message about what coordinates would be valid
    421530            returnval = false;
    422531            allDone = true;
    423         } else if (overlap == PPSTAMP_ERROR) {
     532            foundOverlap = true;
     533            break;
     534        case PPSTAMP_ERROR:
    424535            returnval = false;
    425536            allDone = true;
     537            break;
     538        default:
     539            psError(PS_ERR_PROGRAMMING, false, "findROI returned unexpected value %d\n", overlap);
     540            break;
    426541        }
    427542
    … …  
    430545        if (allDone) {
    431546            view->chip = -1;
    432             view->cell = -1;
    433             break;
    434         } else {
    435             // Remove the transformations from the fpa so that they will be
    436             // recalculated from scratch for the next chip.
    437             // If I don't do this the calculation of the astrometry parameters for the output file
    438             // don't come out correctly. See the code on the false side of the test
    439             // if (fpa->toSky == NULL) in the function pmAstromWCStoFPA
    440             psFree(input->fpa->fromTPA);
    441             psFree(input->fpa->toTPA);
    442             psFree(input->fpa->toSky);
    443             input->fpa->fromTPA = NULL;
    444             input->fpa->toTPA = NULL;
    445             input->fpa->toSky = NULL;
     547            break;
    446548        }
    447549    }
    … …  
    451553    psFree(view);
    452554
     555    if (!foundOverlap) {
     556        fprintf(stderr, "ROI not found in input\n");
     557    }
     558
    453559    return returnval;
    454560}

  • trunk/ppstamp/src/ppstampRegion.c

    r15323 r15363  
    2121    // pmReadoutExtent ignores the col0, row0 of the readout
    2222
    23     int col0 = image->col0;
    24     int row0 = image->row0;
     23    int col0 = 0;
     24    int row0 = 0;
     25   
     26#ifdef notdef
     27    if (image->parent) {
     28        // This is wrong. I'm trying to deal with the strange megacam 'spliced' images
     29        // which appear to have inconsistant CELL.X0 and readout.col0
     30        col0 = image->col0 - image->parent->col0;
     31        row0 = image->row0 - image->parent->row0;
     32    }
     33#endif
    2534
    2635    return psRegionAlloc(col0, col0 + image->numCols,
    … …  
    4453        psFree(roExtent);
    4554    }
    46     // Unlike pmCellExtent we don't add in cell.x0, cell.y0. The appropriate offset
    47     // was included in the roExtent by ppstampReadoutRegion
     55
     56    bool mdok;                          // Status of MD lookup
     57    int cellX0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0"); // Cell x offset
     58    if (!mdok) {
     59        psError(PS_ERR_UNKNOWN, true, "Unable to find CELL.X0.\n");
     60        psFree(cellExtent);
     61        return NULL;
     62    }
     63
     64    int cellY0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0"); // Cell y offset
     65    if (!mdok) {
     66        psError(PS_ERR_UNKNOWN, true, "Unable to find CELL.Y0.\n");
     67        psFree(cellExtent);
     68        return NULL;
     69    }
     70    int xParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
     71    int yParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
     72
     73    if (xParity >= 0) {
     74        cellExtent->x0 += cellX0;
     75        cellExtent->x1 += cellX0;
     76    } else {
     77        float x0 = cellX0 - cellExtent->x1;
     78        float x1 = cellX0 - cellExtent->x0;
     79        cellExtent->x0 = x0;
     80        cellExtent->x1 = x1;
     81    }
     82
     83    if (yParity >= 0) {
     84        cellExtent->y0 += cellY0;
     85        cellExtent->y1 += cellY0;
     86    } else {
     87        float y0 = cellY0 - cellExtent->y1;
     88        float y1 = cellY0 - cellExtent->y0;
     89        cellExtent->y0 = y0;
     90        cellExtent->y1 = y1;
     91    }
    4892
    4993    return cellExtent;

  • trunk/ppstamp/src/ppstampVersion.c

    r15280 r15363  
    2424
    2525
    26 void ppstampVersionMetadata(psMetadata *metadata)
     26void ppstampVersionMetadata(psMetadata *metadata, ppstampOptions *options)
    2727{
    2828    PS_ASSERT_METADATA_NON_NULL(metadata,);
    … …  
    3030    psString pslib = psLibVersionLong();// psLib version
    3131    psString psmodules = psModulesVersionLong(); // psModules version
    32 #ifdef notdef
    33     psString psphot = psphotVersionLong(); // psphot version
    34     psString psastro = psastroVersionLong(); // psastro version
    35     psString ppStats = ppStatsVersionLong(); // ppStats version
    36 #endif
    3732    psString ppstamp = ppstampVersionLong(); // ppstamp version
    3833
    … …  
    4742    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, pslib, "");
    4843    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, psmodules, "");
    49 #ifdef notdef
    50     psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, psphot, "");
    51     psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, psastro, "");
    52     psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, ppStats, "");
    53 #endif
    5444    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, ppstamp, "");
     45    psString roi = NULL;
     46    psStringAppend(&roi, "ppstamp: region of interest: %s", psRegionToString(options->roi));
     47    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY",  PS_META_DUPLICATE_OK, roi, "");
    5548
     49    psFree(roi);
    5650    psFree(head);
    5751    psFree(pslib);
    5852    psFree(psmodules);
    59 #ifdef notdef
    60     psFree(psphot);
    61     psFree(psastro);
    62     psFree(ppStats);
    63 #endif
    6453    psFree(ppstamp);
    6554
Note: See TracChangeset for help on using the changeset viewer.