IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 20271


Ignore:
Timestamp:
Oct 19, 2008, 1:54:21 PM (18 years ago)
Author:
eugene
Message:

save chip-offset values to header

Location:
trunk/psastro/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/psastroFixChips.c

    r16077 r20271  
    66bool psastroFixChips (pmConfig *config, psMetadata *recipe) {
    77
    8   bool status;
    9 
    10   bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS");
    11   if (!status) {
    12       fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
    13   }
    14   if (!fixChips) return true;
    15 
    16   // identify reference astrometry table
    17   // if not defined, correction was not requested; skip step
    18   pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL");
    19   if (!astrom) return true;
    20 
    21   // select the input data sources
    22   pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    23   if (!input) {
    24       psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    25       return false;
    26   }
    27 
    28   // acceptable limits
    29   double pixelTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.TOLERANCE");
    30   if (!status) psAbort ("missing recipe value");
    31 
    32   double angleTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.ANGLE.TOLERANCE");
    33   if (!status) psAbort ("missing recipe value");
    34 
    35   // make sure the astrometry model is loaded.  de-activate all files except PSASTRO.MODEL.
    36   // supply the input concepts so the model is defined using the RA, DEC, POSANGLE of the input
    37   // image.
    38   psFree (astrom->fpa->concepts);
    39   astrom->fpa->concepts = psMemIncrRefCounter (input->fpa->concepts);
    40   pmFPAfileActivate (config->files, false, NULL);
    41   pmFPAfileActivate (config->files, true, "PSASTRO.MODEL");
    42 
    43   pmFPAview *view = pmFPAviewAlloc (0);
    44 
    45   // files associated with the science image
    46   if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    47       psError (PS_ERR_IO, false, "Can't load the astrometry model file");
    48       return false;
    49   }
    50 
    51   // we now have a set of chip solutions and a set of prediction measure the overall
    52   // offset and rotation of the two systems by comparing the chip corners, projected onto
    53   // the focal plane (all 4 to prevent solutions tied to a single corner)
    54 
    55   psVector *xObs = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
    56   psVector *yObs = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
    57   psVector *xRef = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
    58   psVector *yRef = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
    59 
    60   int nPts = 0;
    61 
    62   pmChip *obsChip = NULL;
    63   while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    64       if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
    65 
    66       // set the chip astrometry using the astrom file
    67       pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
    68 
    69       psRegion *region = pmChipPixels (obsChip);
    70       psPlane ptCP, ptFP;
    71 
    72       ptCP.x = region->x0; ptCP.y = region->y0;
    73       psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    74       xObs->data.F32[nPts] = ptFP.x;
    75       yObs->data.F32[nPts] = ptFP.y;
    76       psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    77       xRef->data.F32[nPts] = ptFP.x;
    78       yRef->data.F32[nPts] = ptFP.y;
    79       nPts ++;
    80 
    81       ptCP.x = region->x0; ptCP.y = region->y1;
    82       psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    83       xObs->data.F32[nPts] = ptFP.x;
    84       yObs->data.F32[nPts] = ptFP.y;
    85       psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    86       xRef->data.F32[nPts] = ptFP.x;
    87       yRef->data.F32[nPts] = ptFP.y;
    88       nPts ++;
    89 
    90       ptCP.x = region->x1; ptCP.y = region->y1;
    91       psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    92       xObs->data.F32[nPts] = ptFP.x;
    93       yObs->data.F32[nPts] = ptFP.y;
    94       psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    95       xRef->data.F32[nPts] = ptFP.x;
    96       yRef->data.F32[nPts] = ptFP.y;
    97       nPts ++;
    98 
    99       ptCP.x = region->x1; ptCP.y = region->y0;
    100       psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    101       xObs->data.F32[nPts] = ptFP.x;
    102       yObs->data.F32[nPts] = ptFP.y;
    103       psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    104       xRef->data.F32[nPts] = ptFP.x;
    105       yRef->data.F32[nPts] = ptFP.y;
    106       nPts ++;
    107 
    108       psFree (region);
    109   }
    110   xObs->n = yObs->n = xRef->n = yRef->n = nPts;
     8    bool status;
     9
     10    bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS");
     11    if (!status) {
     12        fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
     13    }
     14    if (!fixChips) return true;
     15
     16    // identify reference astrometry table
     17    // if not defined, correction was not requested; skip step
     18    pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL");
     19    if (!astrom) return true;
     20
     21    // select the input data sources
     22    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
     23    if (!input) {
     24        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     25        return false;
     26    }
     27
     28    // acceptable limits
     29    double pixelTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.TOLERANCE");
     30    if (!status) psAbort ("missing recipe value");
     31
     32    double angleTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.ANGLE.TOLERANCE");
     33    if (!status) psAbort ("missing recipe value");
     34
     35    // make sure the astrometry model is loaded.  de-activate all files except PSASTRO.MODEL.
     36    // supply the input concepts so the model is defined using the RA, DEC, POSANGLE of the input
     37    // image.
     38    psFree (astrom->fpa->concepts);
     39    astrom->fpa->concepts = psMemIncrRefCounter (input->fpa->concepts);
     40    pmFPAfileActivate (config->files, false, NULL);
     41    pmFPAfileActivate (config->files, true, "PSASTRO.MODEL");
     42
     43    pmFPAview *view = pmFPAviewAlloc (0);
     44
     45    // files associated with the science image
     46    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
     47        psError (PS_ERR_IO, false, "Can't load the astrometry model file");
     48        return false;
     49    }
     50
     51    // we now have a set of chip solutions and a set of prediction measure the overall
     52    // offset and rotation of the two systems by comparing the chip corners, projected onto
     53    // the focal plane (all 4 to prevent solutions tied to a single corner)
     54
     55    psVector *xObs = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
     56    psVector *yObs = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
     57    psVector *xRef = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
     58    psVector *yRef = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
     59
     60    int nPts = 0;
     61
     62    pmChip *obsChip = NULL;
     63    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     64        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     65
     66        // set the chip astrometry using the astrom file
     67        pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
     68
     69        psRegion *region = pmChipPixels (obsChip);
     70        psPlane ptCP, ptFP;
     71
     72        ptCP.x = region->x0; ptCP.y = region->y0;
     73        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     74        xObs->data.F32[nPts] = ptFP.x;
     75        yObs->data.F32[nPts] = ptFP.y;
     76        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     77        xRef->data.F32[nPts] = ptFP.x;
     78        yRef->data.F32[nPts] = ptFP.y;
     79        nPts ++;
     80
     81        ptCP.x = region->x0; ptCP.y = region->y1;
     82        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     83        xObs->data.F32[nPts] = ptFP.x;
     84        yObs->data.F32[nPts] = ptFP.y;
     85        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     86        xRef->data.F32[nPts] = ptFP.x;
     87        yRef->data.F32[nPts] = ptFP.y;
     88        nPts ++;
     89
     90        ptCP.x = region->x1; ptCP.y = region->y1;
     91        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     92        xObs->data.F32[nPts] = ptFP.x;
     93        yObs->data.F32[nPts] = ptFP.y;
     94        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     95        xRef->data.F32[nPts] = ptFP.x;
     96        yRef->data.F32[nPts] = ptFP.y;
     97        nPts ++;
     98
     99        ptCP.x = region->x1; ptCP.y = region->y0;
     100        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     101        xObs->data.F32[nPts] = ptFP.x;
     102        yObs->data.F32[nPts] = ptFP.y;
     103        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     104        xRef->data.F32[nPts] = ptFP.x;
     105        yRef->data.F32[nPts] = ptFP.y;
     106        nPts ++;
     107
     108        psFree (region);
     109    }
     110    xObs->n = yObs->n = xRef->n = yRef->n = nPts;
    111111       
    112   psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
     112    psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
    113113 
    114   psVector *mask = psVectorAlloc (nPts, PS_TYPE_U8);
    115   psVectorInit (mask, 0);
    116 
    117   psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    118   stats->clipIter = 1;
    119 
    120   for (int i = 0; i < 3; i++) {
    121       psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef);
    122       psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n);
    123 
    124       psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef);
    125       psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n);
    126   }
    127 
    128   // loop over all chips, select the outliers, and replace the measured astrometry with the model
    129   // the measured transformation above must be applied to make the comparison, and also then applied to the
    130   // model transformation
    131 
    132   psFree (xObs);
    133   psFree (yObs);
    134   psFree (xRef);
    135   psFree (yRef);
    136   psFree (mask);
    137   psFree (stats);
    138 
    139   psFree (view);
    140   view = pmFPAviewAlloc (0);
    141 
    142   while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    143     psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
    144     if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
    145 
    146     // set the chip astrometry using the astrom file
    147     pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
    148 
    149     // bad Astrometry test:  ref pixel or angle outside nominal
    150 
    151     psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
    152     psPlane obsCoord, refCoord, tmpCoord;
    153 
    154     // find location of 0,0 pixel in focal plane coords for this chip
    155     psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel);
    156 
    157     // find location of 0,0 pixel in focal plane coords for ref chip
    158     // apply the global field rotation and offset before comparing
    159     psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel);
    160     psPlaneTransformApply (&refCoord, map, &tmpCoord);
     114    psVector *mask = psVectorAlloc (nPts, PS_TYPE_U8);
     115    psVectorInit (mask, 0);
     116
     117    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     118    stats->clipIter = 1;
     119
     120    for (int i = 0; i < 3; i++) {
     121        psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef);
     122        psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n);
     123
     124        psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef);
     125        psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n);
     126    }
     127
     128    // loop over all chips, select the outliers, and replace the measured astrometry with the model
     129    // the measured transformation above must be applied to make the comparison, and also then applied to the
     130    // model transformation
     131
     132    psFree (xObs);
     133    psFree (yObs);
     134    psFree (xRef);
     135    psFree (yRef);
     136    psFree (mask);
     137    psFree (stats);
     138
     139    psFree (view);
     140    view = pmFPAviewAlloc (0);
     141
     142    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     143        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
     144        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     145
     146        // set the chip astrometry using the astrom file
     147        pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
     148
     149        // bad Astrometry test:  ref pixel or angle outside nominal
     150
     151        psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
     152        psPlane obsCoord, refCoord, tmpCoord;
     153
     154        // find location of 0,0 pixel in focal plane coords for this chip
     155        psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel);
     156
     157        // find location of 0,0 pixel in focal plane coords for ref chip
     158        // apply the global field rotation and offset before comparing
     159        psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel);
     160        psPlaneTransformApply (&refCoord, map, &tmpCoord);
    161161   
    162     psPlane offPixel = {100.0, 0.0, 100.0, 0.0};
    163     psPlane obsOffPt, refOffPt;
    164 
    165     // find location of 0,0 pixel in focal plane coords for this chip
    166     psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel);
    167 
    168     // find location of 0,0 pixel in focal plane coords for ref chip
    169     psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel);
    170     psPlaneTransformApply (&refOffPt, map, &tmpCoord);
     162        psPlane offPixel = {100.0, 0.0, 100.0, 0.0};
     163        psPlane obsOffPt, refOffPt;
     164
     165        // find location of 0,0 pixel in focal plane coords for this chip
     166        psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel);
     167
     168        // find location of 0,0 pixel in focal plane coords for ref chip
     169        psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel);
     170        psPlaneTransformApply (&refOffPt, map, &tmpCoord);
    171171   
    172     double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
    173     double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
    174 
    175     bool badAstrom = false;
    176     badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol;
    177     badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol;
    178     badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
    179 
    180     fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
    181 
    182     if (!badAstrom) continue;
    183 
    184     psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n",
    185               view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
    186 
    187     psFree (obsChip->toFPA);
    188     psFree (obsChip->fromFPA);
    189 
    190     // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
    191     // XXX this only works if toTPA is at most a linear transformation
    192     psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY);
    193     for (int i = 0; i <= refChip->toFPA->x->nX; i++) {
    194         for (int j = 0; j <= refChip->toFPA->x->nY; j++) {
    195             double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
    196             double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
    197             toFPA->x->coeff[i][j] = f1 + f2;
    198 
    199             double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
    200             double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
    201             toFPA->y->coeff[i][j] = g1 + g2;
     172        double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
     173        double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
     174
     175        bool badAstrom = false;
     176        badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol;
     177        badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol;
     178        badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
     179
     180        fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
     181
     182        // XXX for now, just use first readout
     183        pmCell *cell = obsChip->cells->data[0];
     184        pmReadout *readout = cell->readouts->data[0];
     185
     186        // update the header (pull or create local view to entry on readout->analysis)
     187        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     188        if (!updates) {
     189            updates = psMetadataAlloc ();
     190            psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
     191            psFree (updates);
    202192        }
    203     }
    204     toFPA->x->coeff[0][0] += map->x->coeff[0][0];
    205     toFPA->y->coeff[0][0] += map->y->coeff[0][0];
    206 
    207     psRegion *region = pmChipPixels (obsChip);
    208     obsChip->toFPA   = toFPA;
    209     obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50);
    210     psFree (region);
     193
     194        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x);
     195        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y);
     196        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle);
     197
     198        // for successful chips, save the measured offsets in the header
     199        if (!badAstrom) continue;
     200
     201        psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n",
     202                  view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
     203
     204        psFree (obsChip->toFPA);
     205        psFree (obsChip->fromFPA);
     206
     207        // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
     208        // XXX this only works if toTPA is at most a linear transformation
     209        psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY);
     210        for (int i = 0; i <= refChip->toFPA->x->nX; i++) {
     211            for (int j = 0; j <= refChip->toFPA->x->nY; j++) {
     212                double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
     213                double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
     214                toFPA->x->coeff[i][j] = f1 + f2;
     215
     216                double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
     217                double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
     218                toFPA->y->coeff[i][j] = g1 + g2;
     219            }
     220        }
     221        toFPA->x->coeff[0][0] += map->x->coeff[0][0];
     222        toFPA->y->coeff[0][0] += map->y->coeff[0][0];
     223
     224        psRegion *region = pmChipPixels (obsChip);
     225        obsChip->toFPA   = toFPA;
     226        obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50);
     227        psFree (region);
    211228   
    212     // XXX this stuff below should be applied to each readout...
    213     // XXX for now, just use first readout
    214     pmCell *cell = obsChip->cells->data[0];
    215     pmReadout *readout = cell->readouts->data[0];
    216 
    217     // use the new position to re-try the match fit
    218     // select the raw objects for this readout
    219     psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    220     if (rawstars == NULL) { continue; }
    221 
    222     // select the raw objects for this readout
    223     psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    224     if (refstars == NULL) { continue; }
    225 
    226     // the absolute minimum number of stars is 4 (for order = 1)
    227     if ((rawstars->n < 4) || (refstars->n < 4)) {
    228         readout->data_exists = false;
    229         psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)",
    230                   rawstars->n, refstars->n);
    231         continue;
    232     }
    233 
    234     psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars);
    235 
    236     // update the header
    237     psMetadata *updates = psMemIncrRefCounter (psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"));
    238     if (!updates) {
    239         updates = psMetadataAlloc ();
    240     }
    241 
    242     // XXX update the header with info to reflect the failure
    243     if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) {
    244         readout->data_exists = false;
    245         psLogMsg ("psastro", 3, "failed to find a solution\n");
    246         psFree (updates);
    247         continue;
    248     }
    249 
    250     pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL);
    251     psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE | PS_DATA_METADATA, "psastro header stats", updates);
    252     psFree (updates);
    253   }
    254 
    255   psFree (map);
    256   psFree (view);
    257   return true;
     229        // use the new position to re-try the match fit
     230        // select the raw objects for this readout
     231        psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     232        if (rawstars == NULL) { continue; }
     233
     234        // select the raw objects for this readout
     235        psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     236        if (refstars == NULL) { continue; }
     237
     238        // the absolute minimum number of stars is 4 (for order = 1)
     239        if ((rawstars->n < 4) || (refstars->n < 4)) {
     240            readout->data_exists = false;
     241            psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)",
     242                      rawstars->n, refstars->n);
     243            continue;
     244        }
     245
     246        psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars);
     247
     248        // XXX update the header with info to reflect the failure
     249        if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) {
     250            readout->data_exists = false;
     251            psLogMsg ("psastro", 3, "failed to find a solution\n");
     252            continue;
     253        }
     254
     255        pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL);
     256    }
     257
     258    psFree (map);
     259    psFree (view);
     260    return true;
    258261}

  • trunk/psastro/src/psastroFixChipsTest.c

    r16078 r20271  
    201201
    202202        // save WCS and analysis metadata in update header
    203         psMetadata *updates = psMetadataAlloc();
     203        // (pull or create local view to entry on readout->analysis)
     204        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     205        if (!updates) {
     206            updates = psMetadataAlloc ();
     207            psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
     208            psFree (updates);
     209        }
    204210
    205211        // select the raw objects for this readout
    … …  
    228234            readout->data_exists = false;
    229235            psLogMsg ("psastro", 3, "failed to find a solution\n");
    230             psFree (updates);
    231236            continue;
    232237        }
    233238
    234239        pmAstromWriteWCS (updates, fpa, chip, NONLIN_TOL);
    235         psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_DATA_METADATA, "psastro header stats", updates);
    236         psFree (updates);
    237 
    238240    }
    239241
Note: See TracChangeset for help on using the changeset viewer.