IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10830


Ignore:
Timestamp:
Dec 24, 2006, 3:56:53 PM (19 years ago)
Author:
eugene
Message:

finished up the mosaic astrometry code

Location:
trunk/psastro/src
Files:
2 added
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/Makefile.am

    r10716 r10830  
    1212
    1313libpsastro_la_SOURCES = \
    14 psastroArguments.c          \
     14psastroArguments.c          \
     15psastroErrorCodes.c         \
     16psastroVersion.c            \
    1517psastroCleanup.c            \
    1618psastroParseCamera.c        \
    … …  
    1820psastroDataSave.c           \
    1921psastroAstromGuess.c        \
    20 psastroLoadRefstars.c     \
    21 psastroChooseRefstars.c  \
     22psastroLoadRefstars.c       \
     23psastroChooseRefstars.c     \
    2224psastroConvert.c            \
    2325psastroChipAstrom.c         \
    … …  
    2830psastroRefstarSubset.c      \
    2931psastroMosaicAstrom.c       \
    30 psastroMosaicGetGrads.c     \
     32psastroMosaicGradients.c    \
    3133psastroMosaicChipAstrom.c   \
     34psastroMosaicOneChip.c      \
    3235psastroMosaicSetAstrom.c    \
    33 psastroMosaicSetMatch.c     \
    34 psastroMosaicHeaders.c      \
    35 psastroMosaicRescaleChips.c \
    36 psastroErrorCodes.c         \
    37 psastroVersion.c
     36psastroMosaicSetMatch.c
    3837
    3938include_HEADERS = \

  • trunk/psastro/src/psastro.h

    r10716 r10830  
    5656
    5757// mosaic fitting functions
    58 psArray *psastroMosaicGetGrads (pmFPA *fpa, psMetadata *recipe);
     58psArray *psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe);
     59bool psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe);
    5960bool psastroMosaicAstrom (pmConfig *config, psArray *refs);
    60 bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe);
    61 bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe);
     61bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration);
     62bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe, int iteration);
    6263bool psastroMosaicSetAstrom (pmFPA *fpa);
    6364bool psastroMosaicHeaders (pmConfig *config);
    6465bool psastroMosaicRescaleChips (pmFPA *fpa);
    65 
    66 bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale);
    67 psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale);
     66bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration);
    6867
    6968// Return version strings.
    7069psString psastroVersion(void);
    7170psString psastroVersionLong(void);
     71

  • trunk/psastro/src/psastroAstromGuess.c

    r10784 r10830  
    1111    bool newFPA = true;
    1212    bool status = false;
    13     bool isMosaic = false;
    1413    double RAmin = NAN;
    1514    double RAmax = NAN;
    … …  
    2524    DECmin = DECmax = -90;
    2625    RAmin = RAmax = RAminSky = RAmaxSky = 0;
    27 
    28     // is this a mosaic astrometry analysis?
    29     psMetadataLookupStr (&isMosaic, config->arguments, "MOSASTRO");
    3026
    3127    // select the current recipe
    … …  
    4339    }
    4440
    45     double plateScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLATE.SCALE");
    46     if (!status) plateScale = 1.0;
    47     plateScale = 1.0;
     41    // physical pixel scale in microns per pixel
     42    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
     43    if (!status) {
     44        psError(PS_ERR_IO, false, "Failed to lookup pixel scale");
     45        return false;
     46    }
    4847
    4948    pmFPAview *view = pmFPAviewAlloc (0);
    … …  
    5756        pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
    5857
    59         pmAstromReadWCS (fpa, chip, hdu->header, plateScale, isMosaic);
     58        pmAstromReadWCS (fpa, chip, hdu->header, pixelScale);
    6059        if (newFPA) {
    6160            newFPA = false;
    … …  
    6968
    7069        // apply the new WCS guess data to all of the data in the readouts
    71         // XXX should this go into a different function? this would separate WCS interpretation from application
    7270        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
    7371            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    … …  
    8583
    8684                    psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
    87                     psPlaneDistortApply (raw->TP, fpa->toTPA, raw->FP, 0.0, 0.0);
     85                    psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP);
    8886                    psDeproject (raw->sky, raw->TP, fpa->toSky);
    8987
    … …  
    10098                       
    10199                        psProject (tp, raw->sky, fpa->toSky);
    102                         psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0);
     100                        psPlaneTransformApply (fp, fpa->fromTPA, tp);
    103101                        psPlaneTransformApply (ch, chip->fromFPA, fp);
    104102                       

  • trunk/psastro/src/psastroChooseRefstars.c

    r10785 r10830  
    6262
    6363                    psProject (ref->TP, ref->sky, fpa->toSky);
    64                     psPlaneDistortApply (ref->FP, fpa->fromTPA, ref->TP, 0.0, 0.0);
     64                    psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP);
    6565                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
    6666

  • trunk/psastro/src/psastroMosaicAstrom.c

    r10815 r10830  
    11# include "psastro.h"
     2# define NONLIN_TOL 0.001 /* tolerance in pixels */
    23
     4// XXX require this fpa to have multiple chip extensions and a PHU?
    35bool psastroMosaicAstrom (pmConfig *config, psArray *refs) {
    46
    57    bool status;
    6     psArray *grads = NULL;
     8    psArray *gradients = NULL;
    79
    810    // select the current recipe
    … …  
    2022    }
    2123
    22     int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.DISTORTION.ORDER");
    23     if (!status) order = 3.0;
    24 
    2524    pmFPA *fpa = input->fpa;
    2625
    27     psastroMosaicSetMatch (fpa, recipe);
     26    // XXX before we do object matches, we need to fix failed chips
     27    // compare chips with supplied mosaic model
     28    // adjust significant outliers to match model
    2829
    29     // XXX no input distortion model yet; existing fpa distortion is identity: replace it
    30     // XXX make this a test if these are NULL; the load step can NULL them if nothing is loaded
    31     // XXX if we have an input distortion model, the gradient is measuring the error rel to that model
    32     // XXX how does this couple to the individual chip fits?
     30    // given the existing per-chip astrometry, determine matches between raw and ref stars
     31    psastroMosaicSetMatch (fpa, recipe, 0);
     32
     33    // fitted chips will follow the local plate-scale, hiding the distortion
     34    // modify the chip->toFPA scaling to match knowledge about pixel scale
     35    psastroMosaicCommonScale (fpa, recipe);
     36
     37    gradients = psastroMosaicGradients (fpa, recipe);
     38
     39    // allocate mosaic-level polynomial transformation and set masks needed by DVO
     40    int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER");
     41    if (!status) {
     42        psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n");
     43        return false;
     44    }
    3345    psFree (fpa->toTPA);
    34     psFree (fpa->fromTPA);
    35     fpa->toTPA   = psPlaneDistortIdentity (order);
    36     fpa->fromTPA = psPlaneDistortIdentity (order);
     46    fpa->toTPA = psPlaneTransformAlloc (order, order);
     47    for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
     48        for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
     49            if (i + j > order) {
     50                fpa->toTPA->x->mask[i][j] = 1;
     51                fpa->toTPA->y->mask[i][j] = 1;
     52            }
     53        }
     54    }
    3755
    38     grads = psastroMosaicGetGrads (fpa, recipe);
    39 
    40     // fit the measured gradients with the telescope distortion model (3rd order polynomial)
    41     pmAstromFitDistortion (fpa, grads, recipe);
    42     psastroMosaicRescaleChips (fpa);
     56    // fit the measured gradients with the telescope distortion model (polynomial order based on toTPA)
     57    pmAstromFitDistortion (fpa, gradients, recipe);
     58    psFree (gradients);
    4359
    4460    // apply the new distortion terms up and down
    4561    // re-perform the match with a slightly tighter circle
     62    // XXX modify match radius
     63    // XXX set chip.order to 1
    4664    psastroMosaicSetAstrom (fpa);
    47     psastroMosaicSetMatch (fpa, recipe);
     65    psastroMosaicSetMatch (fpa, recipe, 1);
     66    psastroMosaicChipAstrom (fpa, recipe, 1);
    4867
    49     // fit the chips with linear fits
    50     // re-match
    51     psastroMosaicChipAstrom (fpa, recipe);
    52     psastroMosaicSetMatch (fpa, recipe);
     68    // do a second pass on the distortion with improved chip positions and rotations
     69    // XXX do we need to iterate with this step until the distortions don't change?
     70    psastroMosaicCommonScale (fpa, recipe);
     71    gradients = psastroMosaicGradients (fpa, recipe);
     72    pmAstromFitDistortion (fpa, gradients, recipe);
     73    psFree (gradients);
     74   
     75    // now fit the chips under the common distortion with higher-order terms
     76    // XXX modify match radius
     77    // XXX set chip.order to NORDER
     78    psastroMosaicSetAstrom (fpa);
     79    psastroMosaicSetMatch (fpa, recipe, 2);
     80    psastroMosaicChipAstrom (fpa, recipe, 2);
    5381
    54     // re-fit the chips with higher-order fits
    55     // re-match
    56     psastroMosaicChipAstrom (fpa, recipe);
    57     psastroMosaicSetMatch (fpa, recipe);
     82    // save WCS and analysis metadata in update header
     83    // XXX need to add global summary statistics
     84    psMetadata *updates = psMetadataAlloc();
     85    pmAstromWriteBilevelMosaic (updates, fpa, NONLIN_TOL);
     86    psMetadataAdd (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_DATA_METADATA, "psastro header stats", updates);
     87    psFree (updates);
    5888
    59     // re-fit the chips with higher-order fits
    60     psastroMosaicChipAstrom (fpa, recipe);
    61     psastroMosaicHeaders (config);
     89    // update the headers based on the results
     90    // psastroMosaicHeaders (config);
    6291
    6392    return true;
    … …  
    6594
    6695/* coordinate frame hierachy
    67    pixels (on a given readout)
    68    cell
    69    chip
    70    FP (focal plane)
    71    TP (tangent plane)
    72    sky (ra, dec)
    73 */
     96 * pixels (on a given readout)
     97 * cell
     98 * chip
     99 * FP (focal plane)
     100 * TP (tangent plane)
     101 * sky (ra, dec)
     102 */

  • trunk/psastro/src/psastroMosaicChipAstrom.c

    r9574 r10830  
    11# include "psastro.h"
     2# define NONLIN_TOL 0.001 /* tolerance in pixels */
    23
    3 bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe) {
     4bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration) {
    45
    56    pmChip *chip = NULL;
    … …  
    2223                if (! readout->data_exists) { continue; }
    2324
    24                 // select the raw objects for this readout
    25                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    26                 if (rawstars == NULL) { continue; }
     25                // save WCS and analysis metadata in update header
     26                psMetadata *updates = psMetadataAlloc();
    2727
    28                 psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    29                 if (refstars == NULL) { continue; }
     28                psastroMosaicOneChip (chip, readout, recipe, updates, iteration);
    3029
    31                 psArray *match = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH");
    32                 if (match == NULL) { continue; }
    33 
    34                 // need to pass in an update header, sent in from above
    35                 pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, recipe, NULL);
     30                pmAstromWriteBilevelChip (updates, chip, NONLIN_TOL);
     31                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_DATA_METADATA, "psastro header stats", updates);
     32                psFree (updates);
    3633            }
    3734        }
    … …  
    3936    return true;
    4037}
     38
     39/* the iteration value may be 1 or 2.
     40 */

  • trunk/psastro/src/psastroMosaicSetAstrom.c

    r10613 r10830  
    3030       
    3131                    psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
    32                     psPlaneDistortApply (raw->TP, fpa->toTPA, raw->FP, 0.0, 0.0);
     32                    psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP);
    3333                    psDeproject (raw->sky, raw->TP, fpa->toSky);
    3434                }
    … …  
    4141       
    4242                    psProject (ref->TP, ref->sky, fpa->toSky);
    43                     psPlaneDistortApply (ref->FP, fpa->fromTPA, ref->TP, 0.0, 0.0);
     43                    psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP);
    4444                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
    4545                }

  • trunk/psastro/src/psastroMosaicSetMatch.c

    r9574 r10830  
    11# include "psastro.h"
    22
    3 bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe) {
     3bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe, int iteration) {
    44
    55    pmChip *chip = NULL;
    … …  
    77    pmReadout *readout = NULL;
    88    pmFPAview *view = pmFPAviewAlloc (0);
    9 
    10     FILE *f;
    11     char name[128];
     9    char radiusWord[64];
    1210
    1311    FILE *g1 = fopen ("raw.ps.dat", "w");
    1412    FILE *g2 = fopen ("ref.ps.dat", "w");
     13
     14    // use small radius to match stars (assume starting astrometry is good)
     15    bool status = false;
     16    sprintf (radiusWord, "PSASTRO.MOSAIC.RADIUS.N%d", iteration);
     17    double RADIUS = psMetadataLookupF32 (&status, recipe, radiusWord);
     18    if (!status) {
     19        psError(PS_ERR_IO, false, "Failed to lookup matching radius: %s", radiusWord);
     20        return NULL;
     21    }
    1522
    1623    // this loop selects the matched stars for all chips
    … …  
    3744                psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n);
    3845
    39                 // use small radius to match stars (assume starting astrometry is good)
    40                 // XXX should this take a (double radius)?
    41                 psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, recipe);
     46                psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, RADIUS);
    4247                psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n);
    43 
    44                 sprintf (name, "raw.%02d.dat", view->chip);
    45                 f = fopen (name, "w");
    46                 for (int i = 0; i < rawstars->n; i++) {
    47                     pmAstromObj *raw = rawstars->data[i];
    48                     fprintf (f, "%f %f  %f %f  %f %f\n", raw->chip->x, raw->chip->y, raw->FP->x, raw->FP->y, raw->sky->r, raw->sky->d);
    49                 }
    50                 fclose (f);
    51 
    52                 sprintf (name, "ref.%02d.dat", view->chip);
    53                 f = fopen (name, "w");
    54                 for (int i = 0; i < refstars->n; i++) {
    55                     pmAstromObj *ref = refstars->data[i];
    56                     fprintf (f, "%f %f  %f %f  %f %f\n", ref->chip->x, ref->chip->y, ref->FP->x, ref->FP->y, ref->sky->r, ref->sky->d);
    57                 }
    58                 fclose (f);
    5948
    6049                for (int i = 0; i < matches->n; i++) {

  • trunk/psastro/src/psastroOneChip.c

    r10789 r10830  
    2525    psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    2626
     27    // supplied radius is in pixels
     28    double RADIUS = psMetadataLookupF32 (&status, recipe, "PSASTRO.MATCH.RADIUS");
     29    if (!status) {
     30        psError(PS_ERR_IO, false, "Failed to lookup matching radius");
     31        return false;
     32    }
     33
     34    // correct to FP units (physical pixel scale in microns per pixel)
     35    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
     36    if (!status) {
     37        psError(PS_ERR_IO, false, "Failed to lookup pixel scale");
     38        return false;
     39    }
     40    RADIUS *= pixelScale;
     41
    2742    // use small radius to match stars
    28     psArray *match = pmAstromRadiusMatch (rawstars, refstars, recipe);
     43    psArray *match = pmAstromRadiusMatchFP (rawstars, refstars, RADIUS);
    2944    if (stats == NULL) {
    3045        psError(PSASTRO_ERR_UNKNOWN, false, "failed to find radius-matched sources\n");
    … …  
    3247    }
    3348
    34     int nX = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NX");
    35     if (!status) nX = 1;
    36 
    37     int nY = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NY");
    38     if (!status) nY = 1;
    39 
    40     psPlaneTransform *toFPAout = psPlaneTransformAlloc (nX, nY);
    41 
    42     // improved fit for astrometric terms
    43     if (!pmAstromMatchFit (toFPAout, rawstars, refstars, match, recipe, updates)) {
    44         psError(PSASTRO_ERR_DATA, false, "failed to find a valid fitted match solution\n");
     49    // create the output fit model
     50    int order = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.ORDER");
     51    if (!status) {
     52        psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n");
    4553        return false;
    4654    }
    4755    psFree (chip->toFPA);
    48     chip->toFPA = toFPAout;
     56    chip->toFPA = psPlaneTransformAlloc (order, order);
     57    for (int i = 0; i <= chip->toFPA->x->nX; i++) {
     58        for (int j = 0; j <= chip->toFPA->x->nY; j++) {
     59            if (i + j > order) {
     60                chip->toFPA->x->mask[i][j] = 1;
     61                chip->toFPA->y->mask[i][j] = 1;
     62            }
     63        }
     64    }
    4965
    50     // write results
     66    // XXX allow statitic to be set by the user
     67    psStats *fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     68    fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
     69    fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
     70
     71    // improved fit for astrometric terms
     72    pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
     73    if (!results) {
     74        psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n");
     75        return false;
     76    }
     77   
     78    // allowed limits for valid solutions
     79    float maxError = psMetadataLookupF32 (&status, recipe, "PSASTRO.MAX.ERROR");
     80    int minNstar = psMetadataLookupS32 (&status, recipe, "PSASTRO.MIN.NSTAR");
     81
     82    // astError is the average 1D scatter in pixels ('results' are in FPA units = microns)
     83    float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
     84    int astNstar = results->yStats->clippedNvalues;
     85
     86    bool validSolution = true;
     87
     88    // XXX should these result in errors or be handled another way?
     89    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f, Nstars: %d", astError, astNstar);
     90    if (astError > maxError) {
     91        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
     92        validSolution = false;
     93    }
     94    if (astNstar < minNstar) {
     95        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
     96        validSolution = false;
     97    }
     98
     99    // DVO expects NASTRO = 0 if we fail to find a solution
     100    if (validSolution) {
     101        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", astNstar);
     102        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", astError/sqrt(astNstar));
     103    } else {
     104        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", 0);
     105        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", 0.0);
     106    }
     107    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "", astError);
     108    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
     109
     110// write results
    51111    psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    52112
    53113    psFree (match);
    54114    psFree (stats);
     115    psFree (results);
     116    psFree (fitStats);
    55117    psFree (gridStats);
    56118    return true;

  • trunk/psastro/src/psastroUtils.c

    r10613 r10830  
    11# include "psastro.h"
    22# define RENORM 0
     3
     4// fix this to look up the value in the chip concepts
     5static double getChipPixelScale (pmChip *chip) {
     6    return 10.0;
     7}
     8
     9// I have an FPA structure with multiple chips.  we have loaded or measured astrometry for each
     10// chip with independent pixel/degree scaling values.  These will naturally compensate locally
     11// somewhat for the telescope distortion.  to measure the telescope distortion, we need to
     12// force the chips to have the same pixel scale and measure the difference from that solution.
     13// Convert an FPA with disparate pixel scales to a common pixel scale (perhaps depending on the
     14// chip -- eg, TC3)
     15
     16bool psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe) {
     17
     18    // options : use the MIN or MAX chip as global reference or supplied pixel scales
     19    // (microns/pixel), which may depend on the chip
     20
     21    float pixelScaleUse, pixelScale1,  pixelScale2,  pixelScale;
     22
     23    char *option = psMetadataLookupStr (NULL, recipe, "PSASTRO.COMMON.SCALE.OPTION");
     24    if (option == NULL) {
     25        psError(PSASTRO_ERR_DATA, false, "no choice set for common scale option\n");
     26        return false;
     27    }
     28
     29    bool useExternal = true;
     30
     31    // find the min or max scale chip
     32    if (!strcasecmp (option, "MIN") || !strcasecmp (option, "MAX")) {
     33
     34        bool useMax = !strcasecmp (option, "MAX");
     35        pixelScaleUse = (useMax) ? FLT_MIN : FLT_MAX;
     36
     37        for (int i = 0; i < fpa->chips->n; i++) {
     38            pmChip *chip = fpa->chips->data[i];
     39           
     40            pixelScale1 = hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]);
     41            pixelScale2 = hypot (chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
     42            pixelScale = 0.5*(pixelScale1 + pixelScale2);
     43           
     44            pixelScaleUse = (useMax) ? PS_MAX (pixelScale, pixelScaleUse) : PS_MIN (pixelScale, pixelScaleUse);
     45        }
     46        useExternal = false;
     47    }
     48
     49    // rescale each chip by the reference scale
     50    for (int i = 0; i < fpa->chips->n; i++) {
     51        pmChip *chip = fpa->chips->data[i];
     52
     53        psPlaneTransform *toFPA = chip->toFPA;
     54        psPlaneTransform *fromFPA = chip->fromFPA;
     55           
     56        pixelScale1 = hypot (toFPA->x->coeff[1][0], toFPA->x->coeff[0][1]);
     57        pixelScale2 = hypot (toFPA->y->coeff[1][0], toFPA->y->coeff[0][1]);
     58
     59        if (useExternal) {
     60            pixelScaleUse = getChipPixelScale (chip);
     61        }
     62
     63        for (int i = 0; i <= toFPA->x->nX; i++) {
     64            for (int j = 0; j <= toFPA->x->nX; j++) {
     65                toFPA->x->coeff[i][j] *= pixelScaleUse/pixelScale1;
     66                toFPA->y->coeff[i][j] *= pixelScaleUse/pixelScale2;
     67                fromFPA->x->coeff[i][j] *= pixelScale1/pixelScaleUse;
     68                fromFPA->y->coeff[i][j] *= pixelScale2/pixelScaleUse;
     69            }
     70        }
     71
     72    }
     73    return true;
     74}
    375
    476bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars) {
    … …  
    1385
    1486        psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
    15         psPlaneDistortApply (raw->TP, fpa->toTPA, raw->FP, 0.0, 0.0);
     87        psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP);
    1688        psDeproject (raw->sky, raw->TP, fpa->toSky);
    1789    }
Note: See TracChangeset for help on using the changeset viewer.