IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 12523


Ignore:
Timestamp:
Mar 21, 2007, 12:09:22 PM (19 years ago)
Author:
eugene
Message:

fpa field is now determined by the header and/or the chip regions; pixel scale is not needed for readBilevelMosaic, etc; fixed trimsec errors on input/output; fixed errors on map boundaries

Location:
trunk/pswarp/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/pswarp/src/pswarp.h

    r12505 r12523  
    4848pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix);
    4949bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY);
    50 bool pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX, int *nextX);
    51 bool pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY, int *nextY);
     50int pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX);
     51int pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY);
    5252double pswarpMapGridMaxError (pswarpMapGrid *grid);
    5353bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY);

  • trunk/pswarp/src/pswarpDataLoad.c

    r12505 r12523  
    99bool pswarpDataLoad (pmConfig *config) {
    1010
    11     bool status;
    1211    pmChip *chip;
    1312    pmCell *cell;
    1413    pmReadout *readout;
    1514    pmFPAview *view;
    16 
    17     // select the current recipe
    18     psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
    19     if (!recipe) {
    20         psError(PSWARP_ERR_CONFIG, false, "Can't find PSASTRO recipe needed for pixel scale!\n");
    21         return false;
    22     }
    23 
    24     // physical pixel scale in microns per pixel
    25     double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
    26     if (!status) {
    27         psError(PS_ERR_IO, false, "Failed to lookup pixel scale");
    28         return false;
    29     }
    3015
    3116    // select the input data sources
    … …  
    6146    view = pmFPAviewAlloc (0);
    6247
     48    // XXX need to read only the headers for the skycell
     49    // XXX these pmAstromReadBilevel functions seem to be broken
     50
    6351    // find the FPA phu
    64     // XXX wrap the test below into a function
    65     // XXX need to read only the headers for the skycell
    66     // XXX need to optionally load the astrometry datafile
    6752    bool bilevelAstrometry = false;
    6853    pmHDU *phu = pmFPAviewThisPHU (view, astrom->fpa);
    6954    if (phu) {
    70       char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
    71       if (ctype) {
    72         bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    73       }
     55        char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
     56        if (ctype) {
     57            bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
     58        }
    7459    }
    7560    if (bilevelAstrometry) {
    76       pmAstromReadBilevelMosaic (input->fpa, phu->header, pixelScale);
     61        pmAstromReadBilevelMosaic (input->fpa, phu->header);
    7762    }
    7863
    … …  
    8873        pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
    8974        if (bilevelAstrometry) {
    90           pmAstromReadBilevelChip (chip, hdu->header);
     75            pmAstromReadBilevelChip (chip, hdu->header);
    9176        } else {
    92           pmAstromReadWCS (input->fpa, chip, hdu->header, pixelScale);
     77            // we use a default FPA pixel scale of 1.0
     78            pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0);
    9379        }
    9480

  • trunk/pswarp/src/pswarpDefine.c

    r12505 r12523  
    44bool pswarpDefine (pmConfig *config) {
    55
    6     bool status;
    7     psMetadata *recipe = NULL;
    8 
    9     // load the pixel scale from the PSASTRO recipe
    10     recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
    11     if (!recipe) {
    12         psError(PSWARP_ERR_CONFIG, false, "Can't find PSASTRO recipe needed for pixel scale!\n");
    13         return false;
    14     }
    15     // physical pixel scale in microns per pixel
    16     double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
    17     if (!status) {
    18         psError(PS_ERR_IO, false, "Failed to lookup pixel scale");
    19         return false;
    20     }
    21 
    226    // load the PSWARP recipe
    23     recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
     7    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
    248    if (!recipe) {
    259        psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n");
    … …  
    4630    bool bilevelAstrometry = false;
    4731    if (phu) {
    48       char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
    49       if (ctype) {
    50         bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    51       }
     32        char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
     33        if (ctype) {
     34            bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
     35        }
    5236    }
    5337    if (bilevelAstrometry) {
    54       pmAstromReadBilevelMosaic (skycell->fpa, phu->header, pixelScale);
    55       pmAstromReadBilevelChip (chip, hdu->header);
     38        pmAstromReadBilevelMosaic (skycell->fpa, phu->header);
     39        pmAstromReadBilevelChip (chip, hdu->header);
    5640    } else {
    57       pmAstromReadWCS (skycell->fpa, chip, hdu->header, pixelScale);
     41        // we use a default FPA pixel scale of 1.0
     42        pmAstromReadWCS (skycell->fpa, chip, hdu->header, 1.0);
    5843    }
    5944

  • trunk/pswarp/src/pswarpMapGrid.c

    r12505 r12523  
    11# include "pswarp.h"
    22
    3 // construct a grid with superpixel spacing of nXpix, nYpix
    4 // XXX for the moment, ignore readout->cell->chip offsets
     3// pswarpMapGridFromImage builds a set (a grid) of locally-linear maps which convert the source
     4// coordinates (src) to destination coordinates (dest).  we construct a grid with superpixel
     5// spacing of nXpix, nYpix.  The transformation for each grid cell is valid for the superpixel.
     6// The grid over-fills the source image so ever source image pixel is guaranteed to have a map.
    57pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix) {
    68
    … …  
    810    int j, nj;
    911
    10     // XXX I was trying to match the grids too closely
    11     // split the difference of the remainder
    12     // int xMin = 0.5*(src->image->numCols % nXpix);
    13     // int yMin = 0.5*(src->image->numRows % nYpix);
    14     int xMin = 0;
    15     int yMin = 0;
    16 
    17     int nXpts = src->image->numCols / nXpix + 1;
    18     if (src->image->numCols % nXpix) nXpts ++;
    19     int nYpts = src->image->numRows / nYpix + 1;
    20     if (src->image->numRows % nYpix) nYpts ++;
    21 
     12    // start counting from the center of the superpixels
     13    int xMin = 0.5*nXpix;
     14    int yMin = 0.5*nYpix;
     15
     16    // the map is defined for coordinates in the image parent frame.
     17    int Nx = src->image->numCols + src->image->col0;
     18    int Ny = src->image->numRows + src->image->row0;
     19
     20    // allocate an extra superpixel to carry the remainder
     21    int nXpts = Nx / nXpix;
     22    int nYpts = Ny / nYpix;
     23    if (Nx % nXpix) nXpts ++;
     24    if (Ny % nYpix) nYpts ++;
     25
     26    // create the grid of maps
    2227    pswarpMapGrid *grid = pswarpMapGridAlloc (nXpts, nYpts);
    2328
     29    // measure the map for the center of each superpixel
    2430    for (ni = 0, i = xMin; ni < nXpts; i += nXpix, ni++) {
    2531        for (nj = 0, j = yMin; nj < nYpts; j += nYpix, nj++) {
    … …  
    3541}
    3642
     43// set the grid coordinate (gridX,gridY) for the given source image coordinate (ix,iy)
    3744bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY) {
    3845
    39     *gridX = (ix - grid->xMin + 0.5*grid->nXpix) / grid->nXpix;
    40     *gridY = (iy - grid->yMin + 0.5*grid->nYpix) / grid->nYpix;
     46    *gridX = 0.5 + (ix - grid->xMin) / (double) grid->nXpix;
     47    *gridY = 0.5 + (iy - grid->yMin) / (double) grid->nYpix;
    4148    return true;
    4249}
    4350
    44 bool pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX, int *nextX) {
    45 
    46     *nextX = gridX*grid->nXpix + grid->xMin + 0.5*grid->nXpix;
    47     return true;
    48 }
    49 
    50 bool pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY, int *nextY) {
    51 
    52     *nextY = gridY*grid->nYpix + grid->yMin + 0.5*grid->nYpix;
    53     return true;
    54 }
    55 
    56 // measure the max error accumulated in appling one grid point to its neighbors
     51// given the specified grid coordinate (gridX), return the x-coordinate for the source image
     52// corresponding to the next grid cell
     53int pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX) {
     54
     55    int nextX = (gridX + 0.5)*grid->nXpix + grid->xMin;
     56    return nextX;
     57}
     58
     59// given the specified grid coordinate (gridY), return the y-coordinate for the source image
     60// corresponding to the next grid cell
     61int pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY) {
     62
     63    int nextY = (gridY + 0.5)*grid->nYpix + grid->yMin;
     64    return nextY;
     65}
     66
     67// measure the max error accumulated in applying one grid point to its neighbors
     68// XXX double-check this
    5769double pswarpMapGridMaxError (pswarpMapGrid *grid) {
    5870
    … …  
    7688}
    7789
     90// given the source coordinate (inX,inY), return the destination coordinate (outX,outY)
    7891bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY) {
    7992
    … …  
    8497}
    8598
    86 // determine the map for the given pixel from src to dest. pixel is in src coords
     99// determine the (linear) map for the given pixel (ix,iy) from source image (src) to the destination image (dest)
     100// pixel is in src coords. input and output pixel coordinates are in the parent frame of the image (Note that the
     101// astrometric transformations are supplied for the parent image coordinate frame.
    87102bool pswarpMapSetLocalModel (pswarpMap *map, pmReadout *dest, pmReadout *src, int ix, int iy) {
    88103

  • trunk/pswarp/src/pswarpMatchRange.c

    r10958 r12523  
    22
    33# define TEST_MINMAX \
    4     psPlaneTransformApply(FP, chipDest->toFPA, destPix); \
    5     psPlaneTransformApply (TP, fpaDest->toTPA, FP); \
    6     psDeproject (sky, TP, fpaDest->toSky); \
    7     psProject (TP, sky, fpaSrc->toSky); \
    8     psPlaneTransformApply (FP, fpaSrc->fromTPA, TP); \
    9     psPlaneTransformApply (srcPix, chipSrc->fromFPA, FP); \
     4    psPlaneTransformApply(destFP, chipDest->toFPA, destPix); \
     5    psPlaneTransformApply (destTP, fpaDest->toTPA, destFP); \
     6    psDeproject (sky, destTP, fpaDest->toSky); \
     7    psProject (srcTP, sky, fpaSrc->toSky); \
     8    psPlaneTransformApply (srcFP, fpaSrc->fromTPA, srcTP); \
     9    psPlaneTransformApply (srcPix, chipSrc->fromFPA, srcFP); \
     10    fprintf (stderr, "%6.1f,%6.1f -> ", srcPix->x, srcPix->y); \
    1011    *minX = PS_MIN (*minX, srcPix->x); \
    1112    *minY = PS_MIN (*minY, srcPix->y); \
    1213    *maxX = PS_MAX (*maxX, srcPix->x); \
    13     *maxY = PS_MAX (*maxY, srcPix->y);
     14    *maxY = PS_MAX (*maxY, srcPix->y); \
     15    fprintf (stderr, "%4d,%4d - %4d,%4d\n", *minX, *minY, *maxX, *maxY);
    1416
    1517// we are warping from src to dest.  find the max overlapping pixels in the INPUT (src)
    16 // coordinate frame
     18// coordinate frame.  NOTE: these are in the parent pixel frame since the astrometric
     19// transformation refers to the parent frame
    1720bool pswarpMatchRange (int *minX, int *minY, int *maxX, int *maxY, pmReadout *dest, pmReadout *src) {
    1821
    … …  
    3033    pmFPA *fpaDest = chipDest->parent;
    3134
    32     *minX = src->image->numCols;
    33     *minY = src->image->numRows;
    34     *maxX = 0;
    35     *maxY = 0;
     35    *minX = src->image->numCols + src->image->col0 - 1;
     36    *minY = src->image->numRows + src->image->row0 - 1;
     37    *maxX = src->image->col0;
     38    *maxY = src->image->row0;
     39    fprintf (stderr, "%4d,%4d - %4d,%4d\n", *minX, *minY, *maxX, *maxY);
    3640
    3741    // XXX save these as static for speed?
    3842    psPlane *srcPix = psPlaneAlloc();
     43    psPlane *srcFP  = psPlaneAlloc();
     44    psPlane *srcTP  = psPlaneAlloc();
     45
    3946    psPlane *destPix = psPlaneAlloc();
     47    psPlane *destFP  = psPlaneAlloc();
     48    psPlane *destTP  = psPlaneAlloc();
    4049
    41     psPlane *FP = psPlaneAlloc();
    42     psPlane *TP = psPlaneAlloc();
    4350    psSphere *sky = psSphereAlloc();
    4451
    45     destPix->x = 0;
    46     destPix->y = 0;
     52    destPix->x = dest->image->col0;
     53    destPix->y = dest->image->row0;
    4754    TEST_MINMAX;
    4855   
    49     destPix->x = dest->image->numCols;
    50     destPix->y = 0;
     56    destPix->x = dest->image->col0 + dest->image->numCols;
     57    destPix->y = dest->image->row0;
    5158    TEST_MINMAX;
    5259   
    53     destPix->x = dest->image->numCols;
    54     destPix->y = dest->image->numRows;
     60    destPix->x = dest->image->col0 + dest->image->numCols;
     61    destPix->y = dest->image->row0 + dest->image->numRows;
    5562    TEST_MINMAX;
    5663   
    57     destPix->x = 0;
    58     destPix->y = dest->image->numRows;
     64    destPix->x = dest->image->col0;
     65    destPix->y = dest->image->row0 + dest->image->numRows;
    5966    TEST_MINMAX;
    6067   
    61     destPix->x = 0.5*dest->image->numCols;
    62     destPix->y = 0.0;
     68    destPix->x = dest->image->col0 + 0.5*dest->image->numCols;
     69    destPix->y = dest->image->row0;
    6370    TEST_MINMAX;
    6471   
    65     destPix->x = 0.0;
    66     destPix->y = 0.5*dest->image->numRows;
     72    destPix->x = dest->image->col0;
     73    destPix->y = dest->image->row0 + 0.5*dest->image->numRows;
    6774    TEST_MINMAX;
    6875   
    69     destPix->x = 0.5*dest->image->numCols;
    70     destPix->y = dest->image->numRows;
     76    destPix->x = dest->image->col0 + 0.5*dest->image->numCols;
     77    destPix->y = dest->image->row0;
    7178    TEST_MINMAX;
    7279
    73     destPix->x = dest->image->numCols;
    74     destPix->y = 0.5*dest->image->numRows;
     80    destPix->x = dest->image->col0;
     81    destPix->y = dest->image->row0 + 0.5*dest->image->numRows;
    7582    TEST_MINMAX;
    7683
    77     *minX = PS_MAX (*minX, 0);
    78     *minY = PS_MAX (*minY, 0);
    79     *maxX = PS_MIN (*maxX, src->image->numCols);
    80     *maxY = PS_MIN (*maxY, src->image->numRows);
     84    *minX = PS_MAX (*minX, src->image->col0);
     85    *minY = PS_MAX (*minY, src->image->row0);
     86    *maxX = PS_MIN (*maxX, src->image->numCols + src->image->col0);
     87    *maxY = PS_MIN (*maxY, src->image->numRows + src->image->row0);
     88
     89    // demo forward and backward transformation
     90    srcPix->x = *minX;
     91    srcPix->y = *minY;
     92    psPlaneTransformApply(srcFP, chipSrc->toFPA, srcPix);
     93    psPlaneTransformApply (srcTP, fpaSrc->toTPA, srcFP);
     94    psDeproject (sky, srcTP, fpaSrc->toSky);
     95    psProject (destTP, sky, fpaDest->toSky);
     96    psPlaneTransformApply (destFP, fpaDest->fromTPA, destTP);
     97    psPlaneTransformApply (destPix, chipDest->fromFPA, destFP);
     98
     99    fprintf (stderr, "%f,%f -> %f,%f ", srcPix->x, srcPix->y, destPix->x, destPix->y);
     100
     101    psPlaneTransformApply(destFP, chipDest->toFPA, destPix);
     102    psPlaneTransformApply (destTP, fpaDest->toTPA, destFP);
     103    psDeproject (sky, destTP, fpaDest->toSky);
     104    psProject (srcTP, sky, fpaSrc->toSky);
     105    psPlaneTransformApply (srcFP, fpaSrc->fromTPA, srcTP);
     106    psPlaneTransformApply (srcPix, chipSrc->fromFPA, srcFP);
     107
     108    fprintf (stderr, "-> %f,%f ", srcPix->x, srcPix->y);
    81109
    82110    psFree (srcPix);
     111    psFree (srcFP);
     112    psFree (srcTP);
     113
    83114    psFree (destPix);
    84     psFree (FP);
    85     psFree (TP);
     115    psFree (destFP);
     116    psFree (destTP);
     117
    86118    psFree (sky);
    87119

  • trunk/pswarp/src/pswarpTransformReadout_Opt.c

    r12505 r12523  
    55
    66    int minX, minY, maxX, maxY;
    7     int gridXo, gridX, gridY, nextGridX, nextGridY;
     7    int gridXo, gridX, gridY, nextGridXo, nextGridX, nextGridY;
    88    pswarpMap *map = NULL;
    99
    … …  
    1818    psF32 **outData = output->image->data.F32;
    1919
     20    // we need to apply the offset to convert parent coordinates to child coordinates for
     21    // psImagePixelInterpolate below
     22    int inCol0 = input->image->col0;
     23    int inRow0 = input->image->row0;
     24    int outCol0 = output->image->col0;
     25    int outRow0 = output->image->row0;
     26
    2027    // we might want to do the rectangular regions outside of the selection independently
    2128    // psImageInit (output->image, NAN);
    … …  
    2431    pswarpMatchRange (&minX, &minY, &maxX, &maxY, input, output);
    2532
     33    // pswarpMapGridFromImage builds a set of locally-linear maps which convert the
     34    // output coordinates to input coordinates
    2635    pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, 128, 128);
    2736
    … …  
    3140
    3241    pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY);
    33     pswarpMapGridNextGrid_Y (grid, gridY, &nextGridY);
     42    nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
    3443    map = grid->maps[gridX][gridY];
    3544
    36     // Iterate over the output image pixels
     45    assert ((int)(minX - outCol0) >= 0);
     46    assert ((int)(maxX - outCol0) <= output->image->numCols);
     47    assert ((int)(minY - outRow0) >= 0);
     48    assert ((int)(maxY - outRow0) <= output->image->numRows);
     49
     50    gridXo = gridX;
     51    nextGridXo = pswarpMapGridNextGrid_X (grid, gridX);
     52
     53    // Iterate over the output image pixels (parent frame)
    3754    for (int y = minY; y < maxY; y++) {
    3855        if (y >= nextGridY) {
    3956            gridY ++;
    40             pswarpMapGridNextGrid_Y (grid, gridY, &nextGridY);
     57            nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
    4158            map = grid->maps[gridX][gridY];
    4259        }
    4360
    44         gridXo = gridX;
     61        gridX = gridXo;
     62        nextGridX = nextGridXo;
     63        map = grid->maps[gridX][gridY];
    4564        for (int x = minX; x < maxX; x++) {
    4665            if (x >= nextGridX) {
    4766                gridX ++;
    48                 pswarpMapGridNextGrid_X (grid, gridX, &nextGridX);
     67                nextGridX = pswarpMapGridNextGrid_X (grid, gridX);
    4968                map = grid->maps[gridX][gridY];
    5069            }
    … …  
    5372            if (region && region->data.U8[y][x]) continue;
    5473
    55             // XXX double check this 1/2 pixel offset
    56             // XXX subtract 0.5,0.5 from result?
     74            // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
     75            // both are in the parent frames of the input and output images.
    5776            pswarpMapApply (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5);
     77
     78            if (inPix->x - inCol0 < 0) continue;
     79            if (inPix->x - inCol0 >= inImage->numCols) continue;
     80            if (inPix->y - inRow0 < 0) continue;
     81            if (inPix->y - inRow0 >= inImage->numRows) continue;
    5882
    5983            // XXX get interpolation method from the recipe
    6084            // XXX include mask
    6185            // XXX apply scale and offset?
    62             outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
     86            // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates
     87            outData[y-outRow0][x-outCol0] = (psF32)psImagePixelInterpolate(inImage, inPix->x - inCol0, inPix->y - inRow0, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
    6388        }
    64         gridX = gridXo;
    65         pswarpMapGridNextGrid_X (grid, gridX, &nextGridX);
    66         map = grid->maps[gridX][gridY];
    6789    }
    6890
    … …  
    7193    return true;
    7294}
    73 
    74 # if (0)
    75     if (error) {
    76         // Error is actually the variance
    77         outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
    78                                                                                 detector->x,
    79                                                                                 detector->y,
    80                                                                                 mask, 1, NAN);
    81     }
    82     if (error) {
    83         outError->data.F32[y][x] = outError->data.F32[y][x] / SQUARE(scale);
    84     }
    85 # endif
    86 
    87    
Note: See TracChangeset for help on using the changeset viewer.