IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15323


Ignore:
Timestamp:
Oct 16, 2007, 4:00:33 PM (19 years ago)
Author:
bills
Message:

Solved problem with chip and cell boundary calculations by using different
functions than pmChipPixels and pmCellExtents. They don't do their calculations
in image coordinates.

Implemented specifying width/height of the ROI in celestial coordinates.
Removed unneeded WCS calculation methods. (The bug I was chasing was actually
not in ppstamp. It was a bug in ds9

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

Legend:

Unmodified
Added
Removed

  • trunk/ppstamp/src/Makefile.am

    r15280 r15323  
    1717        ppstampCleanup.c \
    1818        ppstampErrorCodes.c \
     19        ppstampMakeStamp.c \
    1920        ppstampOptions.c \
    2021        ppstampParseCamera.c \
    21         ppstampMakeStamp.c \
     22        ppstampRegion.c \
    2223        ppstampVersion.c
    2324

  • trunk/ppstamp/src/ppstamp.h

    r15280 r15323  
    1616#include "ppStats.h"
    1717
    18 #define RECIPE_NAME "PPSTAMP"           // Name of the recipe to use
    1918#define TIMERNAME "ppstamp"             // Name of timer
    2019
    … …  
    2928typedef struct {
    3029    // input arguments
    31     double centerX;
    32     double centerY;
    33     double dX;
    34     double dY;
     30    int    centerX;
     31    int    centerY;
     32    int    dX;
     33    int    dY;
    3534    double centerRA;
    3635    double centerDEC;
    3736    double dRA;
    3837    double dDEC;
    39     bool celestialCenter;       // true if center is in RA/dec
    40     bool celestialRange;        // true if range is in RA/dec
     38    bool   celestialCenter;       // true if center is in RA/dec
     39    bool   celestialRange;        // true if range is in RA/dec
    4140    psString    chipName;
    4241    //
    43     // Calculated Parameters
     42    // Calculated Values
    4443    //
    45     int         chip;
    46     int         cell;
    4744    psRegion    roi;            // roi in chip coordinates
    48     // psPlane     fpaCenter;      // center of ROI in FPA coordinates
    49     // ppstampROI  fpaROI;         // region of interest in FPA coordinates
     45
    5046} ppstampOptions;
    5147
    … …  
    6460bool ppstampMakeStamp(pmConfig *config, ppstampOptions *);
    6561
    66 // Loop over the input
    67 bool ppstampLoop(pmConfig *config, ppstampOptions *options);
    68 
    69 bool ppstampReadout(pmConfig *config, pmFPAview *view);
    70 
    7162// free memory, check for leaks
    7263void ppstampCleanup (pmConfig *config, ppstampOptions *options);
    73 
    74 #ifdef notyet
    75 // calculate stats, including MD5
    76 bool ppstampStats (pmConfig *config, pmChip *chip, const pmFPAview *inputView);
    77 
    78 // write stats to output file
    79 bool ppstampStatsOutput (pmConfig *config);
    80 #endif
    8164
    8265/// Return short version information
    … …  
    9073    );
    9174
    92 
    93 // calculate stats, including MD5
    94 bool ppstampPixelStats (pmConfig *config, const pmFPAview *inputView);
    95 
    96 // calculate stats from headers and concepts
    97 bool ppstampMetadataStats (pmConfig *config);
    98 
    99 void ppstampFileCheck (pmConfig *config);
     75psRegion *ppstampCellRegion(const pmCell *cell);
     76psRegion *ppstampChipRegion(const pmChip *chip);
    10077
    10178#endif

  • trunk/ppstamp/src/ppstampArguments.c

    r15280 r15323  
    44
    55#include "ppstamp.h"
    6 
    7 static void usage (void) {
    8     fprintf(stderr, "USAGE: ppstamp -pixcenter x y    -pixrange dx dy    -chip chipname [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
    9 
    10     fprintf(stderr, "USAGE: ppstamp -pixcenter x y    -celrange dRA dDEC -chip chipname [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
     6#include "ohana.h"
     7
     8static bool get2Angles(int argnum, int *pArgc, char **argv, bool bothDegrees, double *p1, double *p2);
     9static bool get2Ints(int argnum, int *pArgc, char **argv, int *p1, int *p2);
     10
     11// XXX This usage output is kind of busy
     12
     13static void usage (void)
     14{
     15    fprintf(stderr, "\n");
     16
     17    fprintf(stderr, "USAGE: ppstamp -celcenter RA DEC -arcrange dRA dDEC [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
    1118    fprintf(stderr, "       ppstamp -celcenter RA DEC -celrange dRA dDEC [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
    1219    fprintf(stderr, "       ppstamp -celcenter RA DEC -pixrange dx dy    [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n");
    13     fprintf(stderr, "\n");
    14 
    15 #ifdef notyet
     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");
     27    fprintf(stderr, "\n");
    1628    fprintf(stderr, "Optional arguments:\n");
    17     fprintf(stderr, "\t-stats STATS.mdc: Output statistics into STATS.mdc\n");
    18     fprintf(stderr, "\n");
    19     fprintf(stderr, "Input options (single file / file list):\n");
    20     fprintf(stderr, "\n");
    21 #endif
     29    fprintf(stderr, "         -chip chipName    selects chip (only used with -pixcenter)\n");
     30    fprintf(stderr, "\n");
     31
    2232    exit (2);
    23 }
    24 
    25 void get2Doubles(int argnum, int *pArgc, char **argv, double *p1, double *p2) {
    26     if (*pArgc < 3) {
    27         usage();
    28     }
    29     psArgumentRemove(argnum, pArgc, argv);
    30     *p1 = atof(argv[argnum]);
    31     psArgumentRemove(argnum, pArgc, argv);
    32     *p2 = atof(argv[argnum]);
    33     psArgumentRemove(argnum, pArgc, argv);
    3433}
    3534
    … …  
    4847        psString version;
    4948        version = ppstampVersionLong();   fprintf (stdout, "%s\n", version); psFree (version);
    50 // don't need these
    51 //      version = ppStatsVersionLong();   fprintf (stdout, "%s\n", version); psFree (version);
    52 //      version = psphotVersionLong();    fprintf (stdout, "%s\n", version); psFree (version);
    53 //      version = psastroVersionLong();   fprintf (stdout, "%s\n", version); psFree (version);
    5449        version = psModulesVersionLong(); fprintf (stdout, "%s\n", version); psFree (version);
    5550        version = psLibVersionLong();     fprintf (stdout, "%s\n", version); psFree (version);
    … …  
    5853
    5954    // load the site-wide configuration information
    60     pmConfig *config = pmConfigRead(&argc, argv, RECIPE_NAME);
     55    pmConfig *config = pmConfigRead(&argc, argv, NULL);
    6156    if (config == NULL) {
    6257        psErrorStackPrint(stderr, "Can't find site configuration!\n");
    … …  
    6459    }
    6560
    66     (void) pmConfigRecipeOptions (config, RECIPE_NAME);
    67 
    6861    options = ppstampOptionsAlloc();
    6962    *pOptions = options;
    … …  
    7265        gotCenter = true;
    7366        options->celestialCenter = false;
    74         get2Doubles(argnum, &argc, argv, &options->centerX, &options->centerY);
    75     }
     67        psArgumentRemove(argnum, &argc, argv);
     68
     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
    7675    if ((argnum = psArgumentGet(argc, argv, "-celcenter"))) {
    7776        if (gotCenter) {
    … …  
    7978            usage();
    8079        }
    81         get2Doubles(argnum, &argc, argv, &options->centerRA, &options->centerDEC);
    82         options->centerRA  = DEG_TO_RAD(options->centerRA);
    83         options->centerDEC = DEG_TO_RAD(options->centerDEC);
     80        psArgumentRemove(argnum, &argc, argv);
     81
     82        double raDeg, decDeg;
     83        if (!get2Angles(argnum, &argc, argv, false, &raDeg, &decDeg)) {
     84            usage();
     85        }
     86        options->centerRA  = DEG_TO_RAD(raDeg);
     87        options->centerDEC = DEG_TO_RAD(decDeg);
    8488        gotCenter = true;
    8589        options->celestialCenter = true;
    8690    }
     91
    8792    if (!gotCenter) {
    88         usage();
    89     }
     93        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify center of region of interest\n");
     94        usage();
     95    }
     96
    9097    if ((argnum = psArgumentGet(argc, argv, "-pixrange"))) {
    9198        gotRange = true;
    9299        options->celestialRange = false;
    93         get2Doubles(argnum, &argc, argv, &options->dX, &options->dY);
    94     }
    95     if ((argnum = psArgumentGet(argc, argv, "-celrange"))) {
    96         fprintf(stderr, "Specifying range in celestial coordinates is not implemented yet\n");
    97         exit(1);
     100        psArgumentRemove(argnum, &argc, argv);
     101        if (!get2Ints(argnum, &argc, argv, &options->dX, &options->dY)) {
     102            usage();
     103        }
     104    }
     105
     106    if ((argnum = psArgumentGet(argc, argv, "-arcrange"))) {
    98107        if (gotRange) {
    99             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "can't specify both -pixrange and -celrange\n");;
    100             usage();
    101         }
     108            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "specify only one of -pixrange or -arcrange or -celrange\n");;
     109            usage();
     110        }
     111        psArgumentRemove(argnum, &argc, argv);
    102112        gotRange = true;
    103113        options->celestialRange = true;
    104         get2Doubles(argnum, &argc, argv, &options->dRA, &options->dDEC);
    105         options->dRA  = DEG_TO_RAD(options->dRA);
    106         options->dDEC = DEG_TO_RAD(options->dDEC);
    107     }
     114
     115        double deg1, deg2;
     116        if (!get2Angles(argnum, &argc, argv, true, &deg1, &deg2)) {
     117            usage();
     118        }
     119        options->dRA = DEG_TO_RAD(deg1);
     120        options->dDEC = DEG_TO_RAD(deg2);
     121    }
     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    }
     138
    108139    if (!gotRange) {
     140        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify range\n");
    109141        usage();
    110142    }
    … …  
    114146        options->chipName = psStringCopy(argv[argnum]);
    115147        psArgumentRemove(argnum, &argc, argv);
    116     }
    117 
    118     if (!options->celestialCenter && (options->chipName == NULL)) {
    119         fprintf(stderr, "must specify chip name when using -pixcenter\n");
    120         usage();
    121148    }
    122149
    … …  
    125152    if (!status) {
    126153        // TODO: shall we print a specific message?
     154        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify INPUT\n");
    127155        usage ();
    128156    }
    129157
    130     if (argc != 2) usage ();
     158    // finally the output file
     159    if (argc != 2) {
     160        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify OUTPUT\n");
     161        usage ();
     162    }
    131163
    132164    // Add the input and output images (which remain on the command-line) to the arguments list
    … …  
    137169    return config;
    138170}
     171
     172static bool validNumber(char *string)
     173{
     174    char *p = string;
     175    if ((*p == '+') || (*p == '-')) {
     176        p++;
     177    }
     178    return isdigit(*p);
     179}
     180
     181static bool get2Ints(int argnum, int *pArgc, char **argv, int *p1, int *p2)
     182{
     183    if (*pArgc < 2) {
     184        usage();
     185    }
     186
     187    if (!validNumber(argv[argnum])) {
     188        return false;
     189    }
     190
     191    *p1 = atoi(argv[argnum]);
     192    psArgumentRemove(argnum, pArgc, argv);
     193
     194    if (!validNumber(argv[argnum])) {
     195        return false;
     196    }
     197
     198    *p2 = atoi(argv[argnum]);
     199    psArgumentRemove(argnum, pArgc, argv);
     200
     201    return true;
     202}
     203
     204static bool get2Angles(int argnum, int *pArgc, char **argv, bool bothInDegrees, double *p1, double *p2)
     205{
     206    bool rval;
     207
     208    if (*pArgc < 2) {
     209        usage();
     210    }
     211
     212    if (bothInDegrees) {
     213        // both values are angles of arc DD:MM:SS or decimal degrees
     214        rval   = ohana_dms_to_ddd(p1, argv[argnum]);
     215        if (rval) {
     216            rval  = ohana_dms_to_ddd(p1, argv[argnum+1]);
     217        }
     218    } else {
     219        // first value may be in HH:MM:SS
     220        rval = ohana_str_to_radec(p1, p2, argv[argnum], argv[argnum+1]);
     221    }
     222
     223    psArgumentRemove(argnum, pArgc, argv);
     224    psArgumentRemove(argnum, pArgc, argv);
     225
     226    return rval;
     227}

  • trunk/ppstamp/src/ppstampMakeStamp.c

    r15318 r15323  
    1616} ppstampOverlap;
    1717
    18 static bool regionContainsPoint(psRegion *region, psPlane *point);
    19 static bool regionContainsRegion(psRegion *outer, psRegion *inner);
    20 static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, pmCell *inCell,
    21                         psRegion *roi);
    22 
    23 
    24 static ppstampOverlap findCell(pmFPAview *view, ppstampOptions *options, pmFPAfile *input,
    25                                 pmAstromObj *center, pmCell **ppCell)
    26 {
    27     pmCell *cell;
    28 
    29     view->cell = -1;
    30     while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
    31         psRegion *cellExtent = pmCellExtent(cell);
    32         if (regionContainsPoint(cellExtent, center->chip)) {
    33             if (regionContainsRegion(cellExtent, &options->roi)) {
    34                 *ppCell = cell;
    35                 psFree(cellExtent);
    36                 return PPSTAMP_ON;
    37             } else {
    38                 fprintf(stderr, "ROI must be confined to single cell. Partial overlap cell %d\n",
    39                     view->cell);
    40                 psFree(cellExtent);
    41                 return PPSTAMP_PARTIALLY_ON;
    42             }
    43         }
    44         psFree(cellExtent);
    45     }
    46 
    47     // This shouldn't ever happen since ROI is on the chip, it must at least partially overlap
    48     // one of the cells
    49     psError(PS_ERR_PROGRAMMING, false, "findCell couldn't find a cell containing the center\n");
    50 
    51     return PPSTAMP_OFF;
    52 }
     18// convert the input chip's transforms to the output
     19static bool convertWCS(pmHDU *outHDU, pmFPA *outFPA, pmChip *outChip, pmChip *inChip, psRegion *roi)
     20{
     21    pmHDU *hdu = pmHDUFromChip(inChip);
     22    PS_ASSERT_PTR_NON_NULL(hdu, 1)
     23    PS_ASSERT_PTR_NON_NULL(hdu->header, 1)
     24
     25    outChip->toFPA   = psPlaneTransformSetCenter(NULL, inChip->toFPA, roi->x0, roi->y0);
     26    outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50);
     27
     28    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_NONLIN_TOL)) {
     29        psError(PS_ERR_UNKNOWN, false, "Failed to write WCS\n");
     30        return false;
     31    }
     32
     33    return true;
     34}
     35
     36static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, pmCell *inCell, psRegion *roi)
     37{
     38    pmChip    *outChip;
     39    pmFPAview *view = pmFPAviewAlloc(0);
     40
     41    // our output file has a single chip
     42    view->chip = 0;
     43    outChip = pmFPAviewThisChip(view, output->fpa);
     44    psFree(view);
     45
     46    // XXX we probably should copy some concepts from the chip (skipping those that don't apply)
     47    outChip->parent->concepts = psMetadataCopy(outChip->parent->concepts, inChip->parent->concepts);
     48
     49    pmHDU *inHDU  = pmHDUGetHighest(input->fpa, inChip, NULL);
     50    pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);
     51
     52    if (inHDU->header) {
     53        outHDU->header = psMetadataCopy(outHDU->header, inHDU->header);
     54    } else if (!outHDU->header) {
     55        outHDU->header = psMetadataAlloc();
     56    }
     57
     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    }
     71
     72    return true;
     73}
     74
     75// Build the postage stamp output file
    5376
    5477static bool makeStamp(pmConfig *config, ppstampOptions *options, pmFPAfile *input,
    … …  
    82105    while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
    83106        if (!readout->data_exists) {
    84             // XXX if this happens something is wrong
    85107            continue;
    86108        }
    … …  
    92114        }
    93115
    94         // XXX Do we have to do somethign to deal with negative parity
    95116        psImage *subsetImage = psImageSubset(readout->image, options->roi);
    96117        if (subsetImage) {
    … …  
    99120            psFree(subsetImage);
    100121            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.
    101124                // XXX put this into a function in psImage
    102125                // (What I really needed was a function to "create a copy of this portion of an image")
    … …  
    126149
    127150
    128 static bool copyWCS(pmHDU *outHDU, pmFPA *outFPA, pmChip *outChip, pmChip *inChip, psRegion *roi)
    129 {
    130     pmHDU *hdu = pmHDUFromChip(inChip);
    131     PS_ASSERT_PTR_NON_NULL(hdu, 1)
    132     PS_ASSERT_PTR_NON_NULL(hdu->header, 1)
    133 
    134 #ifdef BRUTE_FORCE
    135     outChip->toFPA = inChip->toFPA;
    136     outChip->fromFPA = inChip->fromFPA;
    137     inChip->toFPA = 0;
    138     inChip->fromFPA = 0;
    139 
    140     // This gets us within one pixel of the correct answer
    141     psMetadataItemSupplement(outHDU->header, hdu->header, "CD1_1");
    142     psMetadataItemSupplement(outHDU->header, hdu->header, "CD1_2");
    143     psMetadataItemSupplement(outHDU->header, hdu->header, "CD2_1");
    144     psMetadataItemSupplement(outHDU->header, hdu->header, "CD2_2");
    145     psMetadataItemSupplement(outHDU->header, hdu->header, "CRVAL1");
    146     psMetadataItemSupplement(outHDU->header, hdu->header, "CRVAL2");
    147 
    148     double crpix1 = psMetadataLookupF64(NULL, hdu->header, "CRPIX1");
    149     double crpix2 = psMetadataLookupF64(NULL, hdu->header, "CRPIX2");
    150 
    151     crpix1 -= roi->x0;
    152     crpix2 -= roi->y0;
    153     psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", crpix1);
    154     psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", crpix2);
    155 
    156 
    157 #else
    158 
    159     bool combineTransform = false;
    160     if (combineTransform) {
    161         // translation from output chip coordinates to input chip coordinates
    162         psPlaneTransform *trans = psPlaneTransformAlloc(1, 1);
    163 
    164         trans->x->coeff[0][0] = roi->x0;
    165         trans->x->coeff[0][1] = 0;
    166         trans->x->coeff[1][0] = 1.0;
    167         trans->x->coeff[1][1] = 0;
    168 
    169         trans->y->coeff[0][0] = roi->y0;
    170         trans->y->coeff[0][1] = 1.0;
    171         trans->y->coeff[1][0] = 0;
    172         trans->y->coeff[1][1] = 0;
    173 
    174         // combine the translation with the input's transform
    175         // Note: the region and magic number are not actually used by psPlaneTransformCombine
    176         outChip->toFPA = psPlaneTransformCombine(NULL, trans, inChip->toFPA, *roi, 50);
    177 
    178         // for completeness we set fromFPA, but since we don't use it ...
    179         outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50);
    180 
    181         // psFree(region);
    182         psFree(trans);
     151
     152static bool regionContainsPoint(psRegion *r, psPlane *pt)
     153{
     154    if (pt->x < r->x0)
     155        return false;
     156    if (pt->x >= r->x1)
     157        return false;
     158    if (pt->y < r->y0)
     159        return false;
     160    if (pt->y >= r->y1)
     161        return false;
     162
     163    return true;
     164}
     165
     166// true if the inner region is equal to or completely contained in
     167// the outer region
     168static bool regionContainsRegion(psRegion *outer, psRegion *inner)
     169{
     170    if ((outer->x0 <= inner->x0) &&
     171        (outer->y0 <= inner->y0) &&
     172        (outer->x1 >= inner->x1) &&
     173        (outer->y1 >= inner->y1)) {
     174        return true;
    183175    } else {
    184         outChip->toFPA = psPlaneTransformSetCenter(NULL, inChip->toFPA, roi->x0, roi->y0);
    185         outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50);
    186     }
    187 
    188 
    189     if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_NONLIN_TOL)) {
    190         psError(PS_ERR_UNKNOWN, false, "Failed to write WCS\n");
    191         return false;
    192     }
    193 #endif
    194 
    195     return true;
    196 }
    197 
    198 
    199 static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, pmCell *inCell, psRegion *roi)
    200 {
    201     pmChip    *outChip;
    202     pmFPAview *view = pmFPAviewAlloc(0);
    203 
    204     // our output file has a single chip
    205     view->chip = 0;
    206     outChip = pmFPAviewThisChip(view, output->fpa);
    207     psFree(view);
    208 
    209     outChip->parent->concepts = psMetadataCopy(outChip->parent->concepts, inChip->parent->concepts);
    210 
    211     pmHDU *inHDU  = pmHDUGetHighest(input->fpa, inChip, NULL);
    212     pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);
    213 
    214     if (inHDU->header) {
    215         outHDU->header = psMetadataCopy(outHDU->header, inHDU->header);
    216     } else if (!outHDU->header) {
    217         outHDU->header = psMetadataAlloc();
    218     }
    219 
    220     // steal the input fpa's transforms
    221     output->fpa->toTPA = input->fpa->toTPA;
    222     output->fpa->fromTPA = input->fpa->fromTPA;
    223     output->fpa->toSky = input->fpa->toSky;
    224     // drop references
    225     input->fpa->toTPA = 0;
    226     input->fpa->fromTPA = 0;
    227     input->fpa->toSky = 0;
    228 
    229     if (!copyWCS(outHDU, output->fpa, outChip, inChip, roi)) {
    230         return false;
    231     }
    232 
    233     return true;
    234 }
    235 
    236 
    237 // findROI: find whether the region of interest is completely contained in a cell on a given chip.
    238 //  If so return the cell
    239 ppstampOverlap findROI(ppstampOptions *options, pmFPAview *view, pmFPAfile *input, pmChip *chip,
    240                 pmAstromObj *center, psRegion *roi, pmCell **ppCell)
    241 {
    242     psRegion    *chipExtent = pmChipPixels(chip);
     176        return false;
     177    }
     178}
     179
     180
     181static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
     182{
     183    // convert from sky to FP
     184    psProject(pt->TP, pt->sky, fpa->toSky);
     185    psPlaneTransformApply(pt->FP, fpa->fromTPA, pt->TP);
     186    // convert from FP to chip
     187    psPlaneTransformApply(pt->chip, chip->fromFPA, pt->FP);
     188}
     189
     190static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
     191{
     192    psPlaneTransformApply(pt->FP, chip->toFPA, pt->chip);
     193    psPlaneTransformApply(pt->TP, fpa->toTPA, pt->FP);
     194    psDeproject(pt->sky, pt->TP, fpa->toSky);
     195
     196}
     197
     198static void compareToBox(psRegion *region, psPlane *pt)
     199{
     200    if (pt->x < region->x0)
     201        region->x0 = pt->x;
     202    if (pt->x > region->x1)
     203        region->x1 = pt->x;
     204    if (pt->y < region->y0)
     205        region->y0 = pt->y;
     206    if (pt->y > region->y1)
     207        region->y1 = pt->y;
     208}
     209
     210static void findBoundingBox(ppstampOptions *options, pmFPA *fpa, pmChip *chip, pmAstromObj *center)
     211{
     212    pmAstromObj *pt = pmAstromObjAlloc();
     213
     214    if (!options->celestialCenter) {
     215        // Center was specified in chip coordinates, transform to the sky
     216        chipToSky(center, fpa, chip);
     217    }
     218
     219    // calculate the four corners of the bounding box in sky coordinates, translate them to
     220    // chip coordinates and build the ROI by comparison.
     221
     222    options->roi.x0 = INFINITY;
     223    options->roi.x1 = -INFINITY;
     224    options->roi.y0 = INFINITY;
     225    options->roi.y1 = -INFINITY;
     226   
     227    pt->sky->rErr = 0;
     228    pt->sky->dErr = 0;
     229
     230    pt->sky->r = center->sky->r - options->dRA/2;
     231    pt->sky->d = center->sky->d - options->dDEC/2;
     232    skyToChip(pt, fpa, chip);
     233    compareToBox(&options->roi, pt->chip);
     234
     235    pt->sky->r = center->sky->r + options->dRA/2;
     236    pt->sky->d = center->sky->d - options->dDEC/2;
     237    skyToChip(pt, fpa, chip);
     238    compareToBox(&options->roi, pt->chip);
     239
     240    pt->sky->r = center->sky->r + options->dRA/2;
     241    pt->sky->d = center->sky->d + options->dDEC/2;
     242    skyToChip(pt, fpa, chip);
     243    compareToBox(&options->roi, pt->chip);
     244
     245    pt->sky->r = center->sky->r - options->dRA/2;
     246    pt->sky->d = center->sky->d + options->dDEC/2;
     247    skyToChip(pt, fpa, chip);
     248    compareToBox(&options->roi, pt->chip);
     249
     250    psFree(pt);
     251}
     252
     253
     254// find cell on given chip that contains the region of interest
     255// return the status of the overlap and if the cell completely contains the ROI
     256// set a pointer to reference the cell
     257
     258static ppstampOverlap findCell(pmFPAview *view, ppstampOptions *options, pmFPAfile *input,
     259                                pmAstromObj *center, pmCell **ppCell)
     260{
     261    pmCell *cell;
     262
     263    view->cell = -1;
     264    while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
     265        psRegion *cellExtent = ppstampCellRegion(cell);
     266        if (regionContainsPoint(cellExtent, center->chip)) {
     267            if (regionContainsRegion(cellExtent, &options->roi)) {
     268                *ppCell = cell;
     269                psFree(cellExtent);
     270                return PPSTAMP_ON;
     271            } else {
     272                fprintf(stderr, "ROI must be confined to single cell. Partial overlap cell %d\n",
     273                    view->cell);
     274                fprintf(stderr, "Partial Overlap cell %d\n", view->cell);
     275                fprintf(stderr, "  ROI:         %s\n", psRegionToString(options->roi));
     276                fprintf(stderr, "  Cell Extent: %s\n", psRegionToString(*cellExtent));
     277                psFree(cellExtent);
     278                return PPSTAMP_PARTIALLY_ON;
     279            }
     280        }
     281        psFree(cellExtent);
     282    }
     283
     284    // This shouldn't ever happen since ROI is on the chip, it must at least partially overlap
     285    // one of the cells
     286    psError(PS_ERR_PROGRAMMING, false, "findCell couldn't find a cell containing the center\n");
     287
     288    return PPSTAMP_OFF;
     289}
     290
     291// findROI
     292// 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.
     294
     295static ppstampOverlap findROI(ppstampOptions *options, pmFPAview *view, pmFPAfile *input, pmChip *chip,
     296                       pmAstromObj *center, pmCell **ppCell)
     297{
     298    psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
     299    psRegion    *chipExtent = ppstampChipRegion(chip);
    243300    bool        onChip = false;
    244301    ppstampOverlap   returnval = PPSTAMP_OFF;
    245     psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
    246302
    247303    // set up the astrometry
    … …  
    253309        psError(PS_ERR_UNKNOWN, false, "Failed to Read WCS\n");
    254310        psFree(chipExtent);
    255         psFree(center);
    256311        return PPSTAMP_ERROR;
    257312    }
    … …  
    259314    if (options->celestialCenter) {
    260315
    261         // convert from sky to FP
    262316        center->sky->r = options->centerRA;
    263317        center->sky->d = options->centerDEC;
    264         center->sky->rErr = 0;  // TODO: is this right?
     318        center->sky->rErr = 0;
    265319        center->sky->dErr = 0;
    266320
    267         psProject(center->TP, center->sky, input->fpa->toSky);
    268         psPlaneTransformApply(center->FP, input->fpa->fromTPA, center->TP);
    269 
    270         // convert from FP to chip
     321        skyToChip(center, input->fpa, chip);
    271322        psPlaneTransformApply(center->chip, chip->fromFPA, center->FP);
    272323
    … …  
    277328        }
    278329    } else {
    279         // center specified in pixels, user needs to have specified the chip
    280         if (chipName == NULL) {
    281             psError(PS_ERR_UNKNOWN, true, "Failed to find CHIP.NAME\n");
    282             returnval = PPSTAMP_ERROR;
    283         }
    284         // ppstampArguments insures that options->chipName is not null
    285         if (!strcmp(chipName, options->chipName)) {
     330        // center specified in pixels.
     331        // If the user specified a name of a chip name wait until we get to that one.
     332        // If no chip name was specified, select this one (the first one that had data)
     333        if ((options->chipName == NULL) || !strcmp(chipName, options->chipName)) {
    286334            psLogMsg("ppstampMakeStamp", 3, "Center on chip: %s\n", chipName);
    287             onChip = true;
    288             onChip = true;
    289335            center->chip->x = options->centerX;
    290336            center->chip->y = options->centerY;
     337            center->chip->xErr = 0;
     338            center->chip->yErr = 0;
     339            onChip = true;
    291340        }
    292341    }
    … …  
    294343    if (onChip) {
    295344        if (options->celestialRange) {
    296             // Protect against unimplemented feature
    297             psError(PS_ERR_PROGRAMMING, true, "not ready for Range in celestial coordinates yet.\n");
    298             exit(PS_EXIT_PROG_ERROR);
    299             // find the bounding box in pixel space that encloses the ROI
     345            findBoundingBox(options, input->fpa, chip, center);
    300346        } else {
    301347            int width  = options->dX;
    … …  
    303349
    304350            // calculate the ROI in chip coordinates
    305             roi->x0 = center->chip->x - width / 2;
    306             roi->x1 = roi->x0 + width;
    307             roi->y0 = center->chip->y - height / 2;
    308             roi->y1 = roi->y0 + height;
    309 
    310             if (regionContainsRegion(chipExtent, roi)) {
    311                 returnval = findCell(view, options, input, center, ppCell);
    312             } else {
    313                 fprintf(stderr, "Partial Overlap chip %s\n", chipName);
    314                 fprintf(stderr, "ChipExtent: %s\n", psRegionToString(*chipExtent));
    315                 fprintf(stderr, "ROI:        %s\n", psRegionToString(*roi));
    316 
    317                 returnval = PPSTAMP_PARTIALLY_ON;
    318             }
     351            options->roi.x0 = center->chip->x - width / 2;
     352            options->roi.x1 = options->roi.x0 + width;
     353            options->roi.y0 = center->chip->y - height / 2;
     354            options->roi.y1 = options->roi.y0 + height;
     355        }
     356
     357        if (regionContainsRegion(chipExtent, &options->roi)) {
     358            returnval = findCell(view, options, input, center, ppCell);
     359        } else {
     360            fprintf(stderr, "Partial Overlap chip %s\n", chipName);
     361            fprintf(stderr, "ROI:         %s\n", psRegionToString(options->roi));
     362            fprintf(stderr, "Chip Extent: %s\n", psRegionToString(*chipExtent));
     363
     364            returnval = PPSTAMP_PARTIALLY_ON;
    319365        }
    320366    }
    … …  
    363409
    364410        pmCell *cell;
    365         ppstampOverlap overlap = findROI(options, view, input, chip, center, &options->roi, &cell);
     411        ppstampOverlap overlap = findROI(options, view, input, chip, center, &cell);
    366412
    367413        if (overlap == PPSTAMP_ON) {
    … …  
    409455
    410456
    411 static bool regionContainsPoint(psRegion *r, psPlane *pt)
    412 {
    413     // TODO: should this comparison be done with integers since we
    414     // are dealing with pixels?
    415     // Maybe we should round and truncate before we get here.
    416     if (pt->x < r->x0)
    417         return false;
    418     if (pt->x >= r->x1)
    419         return false;
    420     if (pt->y < r->y0)
    421         return false;
    422     if (pt->y >= r->y1)
    423         return false;
    424 
    425     return true;
    426 }
    427 
    428 // true if the inner region is equal to or completely contained in
    429 // the outer region
    430 static bool regionContainsRegion(psRegion *outer, psRegion *inner)
    431 {
    432     if ((outer->x0 <= inner->x0) &&
    433         (outer->y0 <= inner->y0) &&
    434         (outer->x1 >= inner->x1) &&
    435         (outer->y1 >= inner->y1)) {
    436         return true;
    437     } else {
    438         return false;
    439     }
    440 }

  • trunk/ppstamp/src/ppstampOptions.c

    r15280 r15323  
    1717
    1818    options->celestialCenter = false;
    19     options->centerX         = NAN;
    20     options->centerY         = NAN;
    21     options->centerRA        = NAN;
    22     options->centerDEC       = NAN;
     19    options->centerX         = 0;
     20    options->centerY         = 0;
     21    options->centerRA        = 0;
     22    options->centerDEC       = 0;
    2323    options->celestialRange  = false;
    2424    options->dX   = 0;
Note: See TracChangeset for help on using the changeset viewer.