IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 3, 2010, 8:50:52 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/simtest_nebulous_branches
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/simtest_nebulous_branches

  • branches/simtest_nebulous_branches/pstamp/src/ppstampMakeStamp.c

    r21403 r27840  
    1717
    1818static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);
     19static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);
     20static bool setMaskedToNAN(pmConfig *config, psImage *image, psImage *mask, psImage *variance);
    1921
    2022// convert the input chip's transforms to the output
    21 static bool convertWCS(pmHDU *outHDU, pmFPA *outFPA, pmChip *outChip, pmChip *inChip, psRegion *roi)
    22 {
    23     pmHDU *hdu = pmHDUFromChip(inChip);
    24     PS_ASSERT_PTR_NON_NULL(hdu, 1)
    25     PS_ASSERT_PTR_NON_NULL(hdu->header, 1)
    26 
    27     // Change the reference pixel to account for the change in origin between the stamp and
    28     // the original image
    29     // XXX: shouldn't we be using the center of the stamp instead of the corner?
    30     outChip->toFPA   = psPlaneTransformSetCenter(NULL, inChip->toFPA, (int) roi->x0, (int) roi->y0);
    31 
    32     outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50);
     23static bool convertWCS(pmHDU *outHDU, pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *roi)
     24{
     25    PS_ASSERT_PTR_NON_NULL(outHDU, 1);
     26    PS_ASSERT_PTR_NON_NULL(outHDU->header, 1);
     27
     28    // Center of stamp in chip coordinates
     29    double center_x = (roi->x0 + roi->x1) / 2;
     30    double center_y = (roi->y0 + roi->y1) / 2;
     31    pmAstromObj *check = pmAstromObjAlloc();
     32    check->chip->x = center_x;
     33    check->chip->y = center_y;
     34    chipToSky(check, inFPA, inChip);
     35    double r_before = check->sky->r;
     36    double d_before = check->sky->d;
     37    psLogMsg("ppstampMakeStamp", 2, "Before fit  Center Pixel RA: %f   DEC: %f (degrees)\n",
     38             r_before * PS_DEG_RAD, d_before * PS_DEG_RAD);
    3339
    3440    if (outFPA->toSky->type != PS_PROJ_TAN) {
    35         if (!pmAstromLinearizeTransforms(outFPA, outChip)) {
     41        // we need to make astrometry terms for the output which are centered on the output chip center,
     42        // but we keep the original plate scale
     43        psFree(outFPA->toSky);
     44        outFPA->toSky = psProjectionAlloc (r_before, d_before, inFPA->toSky->Xs, inFPA->toSky->Ys, PS_PROJ_TAN);
     45
     46        if (!pmAstromLinearizeToSky(inFPA, inChip, outFPA, outChip, roi)) {
     47            psFree(check);
    3648            psError(PS_ERR_UNKNOWN, false, "Failed to linearize astrometry\n");
    3749            return false;
    3850        }
    39     }
     51        check->chip->x = (roi->x1 - roi->x0) / 2;
     52        check->chip->y = (roi->y1 - roi->y0) / 2;
     53        chipToSky(check, outFPA, outChip);
     54        double r_after = check->sky->r;
     55        double d_after = check->sky->d;
     56
     57        psLogMsg("ppstampMakeStamp", 2, "After fit:  Center Pixel RA: %f   DEC: %f (degrees)\n",
     58                 r_after * PS_DEG_RAD, d_after * PS_DEG_RAD);
     59        psLogMsg("ppstampMakeStamp", 2, "Error in fit to astrometry %.2f    %.2f (arcseconds)\n",
     60                 (r_after - r_before) *  PS_DEG_RAD * 3600, (d_after - d_before) * PS_DEG_RAD * 3600);
     61
     62        // XXX: should we fail if the fit is bad ??
     63
     64    } else {
     65        outChip->toFPA     = psPlaneTransformSetCenter(NULL, inChip->toFPA, (int) roi->x0, (int) roi->y0);
     66        outChip->fromFPA   = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50);
     67        if (!outChip->fromFPA) {
     68            psError(PS_ERR_UNKNOWN, false, "inversion of toFPA transform failed");
     69            return false;
     70        }
     71    }
     72    psFree(check);
    4073
    4174    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_NONLIN_TOL)) {
    … …  
    4780}
    4881
    49 static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, ppstampOptions *options)
     82static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, ppstampOptions *options, pmAstromObj *center)
    5083{
    5184    pmChip    *outChip;
    … …  
    6093    // since some of the keywords might be duplicated we may not want to copy both
    6194
    62     // copy the fpa concepts
    63     outChip->parent->concepts = psMetadataCopy(outChip->parent->concepts, inChip->parent->concepts);
    64 
    6595    pmHDU *inHDU  = pmHDUFromChip(inChip);
    6696    pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);
    … …  
    72102    }
    73103
     104    // copy the fpa concepts
     105    pmConceptsCopyFPA(output->fpa, input->fpa, false, false);
     106    pmConceptsCopyChip(outChip, inChip, false);
     107    pmCell *outCell = outChip->cells->data[0]; // The only output cell
     108
     109    if (inChip->cells->n == 1) {
     110        pmConceptsCopyCell(outCell, inChip->cells->data[0]);
     111        // Need to fix up the trimsec and biassec to correspond to the output
     112        psMetadataItem *trimsec = psMetadataLookup(outCell->concepts, "CELL.TRIMSEC");
     113        psFree(trimsec->data.V);
     114        trimsec->data.V = NULL;
     115        psMetadataItem *biassec = psMetadataLookup(outCell->concepts, "CELL.BIASSEC");
     116        psFree(biassec->data.V);
     117        biassec->data.V = psListAlloc(NULL);
     118    } else {
     119        psList *inCells = psArrayToList(inChip->cells); // Input cells
     120        pmConceptsAverageCells(outCell, inCells, NULL, NULL, false);
     121        psFree(inCells);
     122    }
    74123
    75124    // If input had WCS convert it for stamp
    76125    if (input->fpa->toSky) {
    77         // steal the input fpa's transforms
    78         output->fpa->toTPA   = input->fpa->toTPA;
    79         output->fpa->fromTPA = input->fpa->fromTPA;
    80         output->fpa->toSky   = input->fpa->toSky;
    81 
    82         // drop the references
    83         input->fpa->toTPA   = 0;
    84         input->fpa->fromTPA = 0;
    85         input->fpa->toSky   = 0;
    86 
    87         if (!convertWCS(outHDU, output->fpa, outChip, inChip, &options->roi)) {
     126        // copy the input fpa's transforms
     127        // XXX: if we change to making multiple stamps per process we'll need to copy
     128        // the transformations
     129        output->fpa->toTPA   = psMemIncrRefCounter(input->fpa->toTPA);
     130        output->fpa->fromTPA = psMemIncrRefCounter(input->fpa->fromTPA);
     131        output->fpa->toSky   = psMemIncrRefCounter(input->fpa->toSky);
     132
     133        if (!convertWCS(outHDU, input->fpa, inChip, output->fpa, outChip, &options->roi)) {
    88134            return false;
    89135        }
    … …  
    92138    }
    93139
     140    psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "RA_DEG", PS_META_REPLACE, "Right Ascension of stamp center", RAD_TO_DEG(center->sky->r));
     141    psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "DEC_DEG", PS_META_REPLACE, "Declination of stamp center", RAD_TO_DEG(center->sky->d));
     142
    94143    ppstampVersionMetadata(outHDU->header, options);
    95144
    … …  
    99148// Extract pixels from an image. If part of the bounds of the region are
    100149// partially outside of the input those pixels are set to zero.
    101 psImage *extractStampOld(psImage *image, psRegion region)
    102 {
    103     int width  = region.x1 - region.x0;
    104     int height = region.y1 - region.y0;
     150static psImage *extractStamp(psImage *image, psRegion region, double value)
     151{
     152    int width  = region.x1 - region.x0 + 0.5;
     153    int height = region.y1 - region.y0 + 0.5;
    105154
    106155    if (width < 0) {
     156        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "negative width\n");
    107157        return NULL;
    108158    }
    109159    if (height < 0) {
     160        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "negative height\n");
    110161        return NULL;
    111162    }
    … …  
    131182
    132183    int copyWidth  = lastX - srcX;
    133     int rightBlank = width - copyWidth - leftBlank;
    134     if (copyWidth <= 0) {
    135         return NULL;
    136     }
    137 
    138     psImage *output = psImageAlloc(width, height, image->type.type);
    139 
    140     psElemType  type = image->type.type;
    141     psS32       elementSize =  PSELEMTYPE_SIZEOF(type);
    142     int         copySize = copyWidth * elementSize;
    143 
    144     for (int dstY=0 ; dstY < height; srcY++, dstY++) {
    145         psU8 *pdst = output->data.U8[dstY];
    146 
    147         if ((srcY < image->row0) || (srcY >= lastY)) {
    148             // zero the row
    149             memset(pdst, 0, width*elementSize);
    150         } else {
    151             psU8 *psrc = image->data.U8[srcY]  + (srcX * elementSize);
    152 
    153             if (leftBlank) {
    154                 memset(pdst, 0, leftBlank*elementSize);
    155             }
    156 
    157             // copy copyWidth pixels from srcX,srcY to dstX, dstY
    158             pdst += dstX*elementSize;
    159             memcpy(pdst, psrc, copySize);
    160 
    161             if (rightBlank) {
    162                 pdst += copyWidth * elementSize;
    163                 memset(pdst, 0, rightBlank);
    164             }
    165         }
    166     }
    167 
    168 
    169     return output;
    170 }
    171 // Extract pixels from an image. If part of the bounds of the region are
    172 // partially outside of the input those pixels are set to zero.
    173 static psImage *extractStamp(psImage *image, psRegion region, double value)
    174 {
    175     int width  = region.x1 - region.x0;
    176     int height = region.y1 - region.y0;
    177 
    178     if (width < 0) {
    179         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "negative width\n");
    180         return NULL;
    181     }
    182     if (height < 0) {
    183         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "negative height\n");
    184         return NULL;
    185     }
    186 
    187     int srcX  = region.x0;
    188     int srcY  = region.y0;
    189     int lastX = region.x1;
    190     int lastY = region.y1;
    191     if (lastX > (image->col0 + image->numCols)) {
    192         lastX = image->col0 + image->numCols;
    193     }
    194     if (lastY > (image->row0 + image->numRows)) {
    195         lastY = image->row0 + image->numRows;
    196     }
    197 
    198     int leftBlank = 0;
    199     int dstX  = 0;
    200     if (srcX < image->col0) {
    201         leftBlank = image->col0 - srcX;
    202         dstX += leftBlank;
    203         srcX = 0;
    204     }
    205 
    206     int copyWidth  = lastX - srcX;
    207 //    int rightBlank = width - copyWidth - leftBlank;
    208184    if (copyWidth <= 0) {
    209185        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "copyWidth is invalid: %d\n", copyWidth);
    210186        return NULL;
     187    }
     188    if (copyWidth > width) {
     189        copyWidth = width;
    211190    }
    212191
    … …  
    230209    if (skip > 0) {
    231210        dstY = skip;
    232         height -= skip;
    233211        srcY = image->row0;
    234212    }
    … …  
    237215            break;
    238216        }
    239         psU8 *pdst = output->data.U8[dstY];
     217        // copy copyWidth pixels from srcX,srcY to dstX, dstY
     218        psU8 *pdst = output->data.U8[dstY] + (dstX * elementSize);
    240219
    241220        psU8 *psrc = image->data.U8[srcY]  + (srcX * elementSize);
    242221
    243         // copy copyWidth pixels from srcX,srcY to dstX, dstY
    244         pdst += dstX*elementSize;
    245222        memcpy(pdst, psrc, copySize);
    246223    }
    … …  
    251228// Build the postage stamp output file
    252229
    253 static bool makeStamp(pmConfig *config, ppstampOptions *options, pmFPAfile *input,
    254                 pmChip *inChip, pmFPAview *view)
     230static int makeStamp(pmConfig *config, ppstampOptions *options, pmFPAfile *input,
     231                pmChip *inChip, pmFPAview *view, pmAstromObj *center)
    255232{
    256233    int status = false;
    … …  
    259236    if (!output) {
    260237        psError(PS_ERR_UNKNOWN, false, "Can't find output data\n");
    261         return false;
     238        return PS_EXIT_DATA_ERROR;
    262239    }
    263240    char *fpaName = psMetadataLookupStr(NULL, input->fpa->concepts, "FPA.OBS"); // Name of FPA
    … …  
    283260        pmFPAfile *mosaic = ppstampBuildMosaic(config, input, view);
    284261        if (mosaic == NULL) {
    285             return false;
     262            return PS_EXIT_UNKNOWN_ERROR;
    286263        }
    287264        srcFile = mosaic;
    … …  
    319296        }
    320297
    321         outReadout->image = extractStamp(readout->image, extractRegion,  0);
     298        outReadout->image = extractStamp(readout->image, extractRegion,  NAN);
    322299        if (!outReadout->image) {
    323300            psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp image\n");
    324301            status = false;
    325302            break;
     303        }
     304        if (readout->variance) {
     305            outReadout->variance = extractStamp(readout->variance, extractRegion,  NAN);
     306            if (!outReadout->variance) {
     307                psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp weight image\n");
     308                status = false;
     309                break;
     310            }
    326311        }
    327312        if (readout->mask) {
    … …  
    332317                break;
    333318            }
    334         }
    335         if (readout->variance) {
    336             outReadout->variance = extractStamp(readout->variance, extractRegion,  0);
    337             if (!outReadout->variance) {
    338                 psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp weight image\n");
    339                 status = false;
    340                 break;
     319
     320            if (!setMaskedToNAN(config, outReadout->image, outReadout->mask, outReadout->variance)) {
     321                 psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp mask image\n");
     322                 status = false;
     323                 break;
    341324            }
    342325        }
    … …  
    353336
    354337    if (status) {
    355         status = copyMetadata(output, input, inChip, options);
    356     }
    357     return status;
     338        status = copyMetadata(output, input, inChip, options, center);
     339    }
     340    return status ? PS_EXIT_SUCCESS : PS_EXIT_UNKNOWN_ERROR;
    358341}
    359342
    … …  
    389372
    390373
    391 static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
    392 {
    393     // convert from sky to TP to FP
    394     psProject(pt->TP, pt->sky, fpa->toSky);
     374static void TPToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
     375{
    395376    psPlaneTransformApply(pt->FP, fpa->fromTPA, pt->TP);
    396377    // convert from FP to chip
    397378    psPlaneTransformApply(pt->chip, chip->fromFPA, pt->FP);
     379}
     380static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
     381{
     382    // convert from sky to TP to FP
     383    psProject(pt->TP, pt->sky, fpa->toSky);
     384    TPToChip(pt, fpa, chip);
    398385}
    399386
    … …  
    405392    psPlaneTransformApply(pt->TP, fpa->toTPA, pt->FP);
    406393    psDeproject(pt->sky, pt->TP, fpa->toSky);
    407 
    408394}
    409395
    … …  
    426412    pmAstromObj *pt = pmAstromObjAlloc();
    427413
    428     if (!options->roip.celestialCenter) {
    429         // Center was specified in chip coordinates, transform to the sky
    430         chipToSky(center, fpa, chip);
    431     }
    432 
    433414    // calculate the four corners of the bounding box in sky coordinates, translate them to
    434415    // chip coordinates and build the ROI by comparison.
    … …  
    439420    options->roi.y1 = -INFINITY;
    440421
    441     pt->sky->rErr = 0;
    442     pt->sky->dErr = 0;
    443 
    444     pt->sky->r = center->sky->r - options->roip.dRA/2;
    445     pt->sky->d = center->sky->d - options->roip.dDEC/2;
    446     skyToChip(pt, fpa, chip);
     422    double dx = 0.5 * options->roip.dRA  / fpa->toSky->Xs;
     423    double dy = 0.5 * options->roip.dDEC / fpa->toSky->Ys;
     424
     425    if (dx > 5000) {
     426        fprintf(stderr, "requested width %f too large reducing to 5000\n", dx);
     427        dx = 5000;
     428    }
     429    if (dy > 5000) {
     430        fprintf(stderr, "requested height %f too large reducing to 5000\n", dy);
     431        dy = 5000;
     432    }
     433
     434    pt->TP->xErr = 0;
     435    pt->TP->yErr = 0;
     436
     437    pt->TP->x = center->TP->x - dx;
     438    pt->TP->y = center->TP->y - dy;
     439    TPToChip(pt, fpa, chip);
    447440    compareToBox(&options->roi, pt->chip);
    448441
    449     pt->sky->r = center->sky->r + options->roip.dRA/2;
    450     pt->sky->d = center->sky->d - options->roip.dDEC/2;
    451     skyToChip(pt, fpa, chip);
     442    pt->TP->x = center->TP->x + dx;
     443    pt->TP->y = center->TP->y - dy;
     444    TPToChip(pt, fpa, chip);
    452445    compareToBox(&options->roi, pt->chip);
    453446
    454     pt->sky->r = center->sky->r + options->roip.dRA/2;
    455     pt->sky->d = center->sky->d + options->roip.dDEC/2;
    456     skyToChip(pt, fpa, chip);
     447    pt->TP->x = center->TP->x + dx;
     448    pt->TP->y = center->TP->y + dy;
     449    TPToChip(pt, fpa, chip);
    457450    compareToBox(&options->roi, pt->chip);
    458451
    459     pt->sky->r = center->sky->r - options->roip.dRA/2;
    460     pt->sky->d = center->sky->d + options->roip.dDEC/2;
    461     skyToChip(pt, fpa, chip);
     452    pt->TP->x = center->TP->x - dx;
     453    pt->TP->y = center->TP->y + dy;
     454    TPToChip(pt, fpa, chip);
    462455    compareToBox(&options->roi, pt->chip);
    463456
    … …  
    532525            center->chip->yErr = 0;
    533526            onChip = true;
     527            chipToSky(center, input->fpa, chip);
    534528        }
    535529    }
    … …  
    541535            int width  = options->roip.dX;
    542536            int height = options->roip.dY;
     537            if (width > 5000) {
     538                fprintf(stderr, "requested width %d too large reducing to 5000\n", width);
     539                width = 5000;
     540            }
     541            if (height > 5000) {
     542                fprintf(stderr, "requested height %d too large reducing to 5000\n", height);
     543                height = 5000;
     544            }
    543545
    544546            // calculate the ROI in chip coordinates
    … …  
    567569}
    568570
    569 bool ppstampMakeStamp (pmConfig *config, ppstampOptions *options)
     571int ppstampMakeStamp (pmConfig *config, ppstampOptions *options)
    570572{
    571573    bool        status = false;
    572     bool        returnval = false;;
     574    int        returnval = PS_EXIT_SUCCESS;;
    573575    bool        foundOverlap = false;
    574576
    … …  
    576578    if (!status) {
    577579        psError(PS_ERR_UNKNOWN, true, "Can't find input file!\n");
    578         return false;
     580        return PS_EXIT_DATA_ERROR;
    579581    }
    580582
    … …  
    584586    } else if (astrom->camera != input->camera) {
    585587        psError(PS_ERR_UNKNOWN, true, "Input camera and astrometry camera do not match");
    586         return false;
     588        return PS_EXIT_CONFIG_ERROR;
    587589    }
    588590
    … …  
    593595        psError(PS_ERR_UNKNOWN, false, "Failed to load input.");
    594596        psFree (view);
    595         return false;
     597        return PS_EXIT_DATA_ERROR;
    596598    }
    597599    bool bilevelAstrometry  = false;
    … …  
    607609            psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel mosaic astrometry for input FPA.");
    608610            psFree(view);
    609             return false;
     611            return PS_EXIT_DATA_ERROR;
    610612        }
    611613    }
    … …  
    635637        case PPSTAMP_ON:
    636638        case PPSTAMP_PARTIALLY_ON:
    637             returnval = makeStamp(config, options, input, chip, view);
     639            returnval = makeStamp(config, options, input, chip, view, center);
    638640            allDone = true;
    639641            foundOverlap = true;
    640642            break;
    641643        case PSTAMP_ERROR:
    642             returnval = false;
     644            returnval = PS_EXIT_UNKNOWN_ERROR;
    643645            allDone = true;
    644646            break;
    … …  
    660662    psFree(view);
    661663
    662     if (!foundOverlap) {
     664    if (!foundOverlap && (returnval == PS_EXIT_SUCCESS)) {
    663665        fprintf(stderr, "ROI not found in input\n");
     666        returnval = PSTAMP_NO_OVERLAP;
    664667    }
    665668
    … …  
    668671
    669672
     673static bool setMaskedToNAN(pmConfig *config, psImage *image, psImage *mask, psImage *variance)
     674{
     675    bool status;
     676    psMetadata *masks = psMetadataLookupMetadata(&status, config->recipes, "MASKS");
     677    if (!status) {
     678        psError(PM_ERR_CONFIG, false, "failed to lookup MASKS in recipes\n");
     679        return false;
     680    }
     681    // we set anything masked to NAN except if CONV.POOR is the only bit set
     682    // First check the old value
     683    psU32 convPoor = psMetadataLookupU32(&status, masks, "POOR.WARP");
     684    if (!status) {
     685        convPoor = psMetadataLookupU32(&status, masks, "CONV.POOR");
     686        if (!status) {
     687            psError(PM_ERR_CONFIG, false, "failed to lookup mask value for CONV.POOR in recipes\n");
     688            return false;
     689        }
     690    }
     691    psU32 maskMask = ~convPoor;
     692
     693    double exciseValue;
     694    if (image->type.type == PS_TYPE_U16) {
     695        exciseValue = 0xffff;
     696    } else if (image->type.type == PS_TYPE_F32) {
     697        exciseValue = NAN;
     698    } else {
     699         psError(PS_ERR_PROGRAMMING, true, "unexpected image type: %d\n", image->type.type);
     700        return false;
     701    }
     702    long numExcised = 0;
     703    for (int y=0; y<image->numRows; y++) {
     704        for (int x=0; x<image->numCols; x++) {
     705            psU16 maskVal = psImageGet(mask, x, y);
     706            if (maskVal & maskMask) {
     707                numExcised++;
     708                psImageSet(image, x, y, exciseValue);
     709                if (variance) {
     710                    psImageSet(variance, x, y, exciseValue);
     711                }
     712            }
     713        }
     714    }
     715    fprintf(stderr, "excised %ld masked pixels\n", numExcised);
     716
     717    return true;
     718}
Note: See TracChangeset for help on using the changeset viewer.