Changeset 16077

Timestamp:

Jan 14, 2008, 4:56:48 PM (18 years ago)

Author:

eugene

Message:

substantial rework: calculate the average rotation and offset, apply to the TPA

File:

• trunk/psastro/src/psastroFixChips.c (modified) (5 diffs)

trunk/psastro/src/psastroFixChips.c

@@ -49 +49 @@
   }
 
-  // loop over all chips, replace input astrometry elements with those from astrom 
+  // we now have a set of chip solutions and a set of prediction measure the overall
+  // offset and rotation of the two systems by comparing the chip corners, projected onto
+  // the focal plane (all 4 to prevent solutions tied to a single corner)
+
+  psVector *xObs = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
+  psVector *yObs = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
+  psVector *xRef = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
+  psVector *yRef = psVectorAllocEmpty (4*input->fpa->chips->n, PS_TYPE_F32);
+
+  int nPts = 0;
+
   pmChip *obsChip = NULL;
+  while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
+      if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
+
+      // set the chip astrometry using the astrom file
+      pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
+
+      psRegion *region = pmChipPixels (obsChip);
+      psPlane ptCP, ptFP;
+
+      ptCP.x = region->x0; ptCP.y = region->y0;
+      psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
+      xObs->data.F32[nPts] = ptFP.x;
+      yObs->data.F32[nPts] = ptFP.y;
+      psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
+      xRef->data.F32[nPts] = ptFP.x;
+      yRef->data.F32[nPts] = ptFP.y;
+      nPts ++;
+
+      ptCP.x = region->x0; ptCP.y = region->y1;
+      psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
+      xObs->data.F32[nPts] = ptFP.x;
+      yObs->data.F32[nPts] = ptFP.y;
+      psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
+      xRef->data.F32[nPts] = ptFP.x;
+      yRef->data.F32[nPts] = ptFP.y;
+      nPts ++;
+
+      ptCP.x = region->x1; ptCP.y = region->y1;
+      psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
+      xObs->data.F32[nPts] = ptFP.x;
+      yObs->data.F32[nPts] = ptFP.y;
+      psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
+      xRef->data.F32[nPts] = ptFP.x;
+      yRef->data.F32[nPts] = ptFP.y;
+      nPts ++;
+
+      ptCP.x = region->x1; ptCP.y = region->y0;
+      psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
+      xObs->data.F32[nPts] = ptFP.x;
+      yObs->data.F32[nPts] = ptFP.y;
+      psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
+      xRef->data.F32[nPts] = ptFP.x;
+      yRef->data.F32[nPts] = ptFP.y;
+      nPts ++;
+
+      psFree (region);
+  }
+  xObs->n = yObs->n = xRef->n = yRef->n = nPts;
+        
+  psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
+  
+  psVector *mask = psVectorAlloc (nPts, PS_TYPE_U8);
+  psVectorInit (mask, 0);
+
+  psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+  stats->clipIter = 1;
+
+  for (int i = 0; i < 3; i++) {
+      psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef);
+      psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n);
+
+      psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef);
+      psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n);
+  }
+
+  // loop over all chips, select the outliers, and replace the measured astrometry with the model
+  // the measured transformation above must be applied to make the comparison, and also then applied to the 
+  // model transformation
+
+  psFree (xObs);
+  psFree (yObs);
+  psFree (xRef);
+  psFree (yRef);
+  psFree (mask);
+  psFree (stats);
+
+  psFree (view);
+  view = pmFPAviewAlloc (0);
+
   while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
@@ -61 +150 @@
 
     psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
-    psPlane obsCoord, refCoord;
+    psPlane obsCoord, refCoord, tmpCoord;
 
     // find location of 0,0 pixel in focal plane coords for this chip
@@ -67 +156 @@
 
     // find location of 0,0 pixel in focal plane coords for ref chip
-    psPlaneTransformApply (&refCoord, refChip->toFPA, &refPixel);
+    // apply the global field rotation and offset before comparing
+    psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel);
+    psPlaneTransformApply (&refCoord, map, &tmpCoord);
     
-    psPlane offPixel = {0.0, 0.0, 0.0, 0.0};
+    psPlane offPixel = {100.0, 0.0, 100.0, 0.0};
     psPlane obsOffPt, refOffPt;
 
@@ -76 +167 @@
 
     // find location of 0,0 pixel in focal plane coords for ref chip
-    psPlaneTransformApply (&refOffPt, refChip->toFPA, &offPixel);
+    psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel);
+    psPlaneTransformApply (&refOffPt, map, &tmpCoord);
     
-    double obsAngle = atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
-    double refAngle = atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
+    double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
+    double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
 
     bool badAstrom = false;
@@ -86 +178 @@
     badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
 
+    fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
+
     if (!badAstrom) continue;
+
+    psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n",
+              view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
 
     psFree (obsChip->toFPA);
     psFree (obsChip->fromFPA);
 
-    // supply astrometry from model
-    obsChip->toFPA   = psMemIncrRefCounter (refChip->toFPA);
-    obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA);
-  }
-
-  psFree (input->fpa->toSky);
-  psFree (input->fpa->toTPA);
-  psFree (input->fpa->fromTPA);
-  input->fpa->toSky   = psMemIncrRefCounter (astrom->fpa->toSky);
-  input->fpa->toTPA   = psMemIncrRefCounter (astrom->fpa->toTPA);
-  input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA);
-
-  psMetadata *updates = psMetadataAlloc();
-  pmAstromWriteBilevelMosaic (updates, input->fpa, NONLIN_TOL);
-  psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
-  psFree (updates);
-
+    // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
+    // XXX this only works if toTPA is at most a linear transformation
+    psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY);
+    for (int i = 0; i <= refChip->toFPA->x->nX; i++) {
+        for (int j = 0; j <= refChip->toFPA->x->nY; j++) {
+            double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
+            double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
+            toFPA->x->coeff[i][j] = f1 + f2;
+
+            double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
+            double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
+            toFPA->y->coeff[i][j] = g1 + g2;
+        }
+    }
+    toFPA->x->coeff[0][0] += map->x->coeff[0][0];
+    toFPA->y->coeff[0][0] += map->y->coeff[0][0];
+
+    psRegion *region = pmChipPixels (obsChip);
+    obsChip->toFPA   = toFPA;
+    obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50);
+    psFree (region);
+    
+    // XXX this stuff below should be applied to each readout...
+    // XXX for now, just use first readout
+    pmCell *cell = obsChip->cells->data[0];
+    pmReadout *readout = cell->readouts->data[0];
+
+    // use the new position to re-try the match fit
+    // select the raw objects for this readout
+    psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
+    if (rawstars == NULL) { continue; }
+
+    // select the raw objects for this readout
+    psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
+    if (refstars == NULL) { continue; }
+
+    // the absolute minimum number of stars is 4 (for order = 1)
+    if ((rawstars->n < 4) || (refstars->n < 4)) {
+        readout->data_exists = false;
+        psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", 
+                  rawstars->n, refstars->n);
+        continue;
+    } 
+
+    psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars);
+
+    // update the header
+    psMetadata *updates = psMemIncrRefCounter (psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"));
+    if (!updates) {
+        updates = psMetadataAlloc ();
+    }
+
+    // XXX update the header with info to reflect the failure
+    if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) {
+        readout->data_exists = false;
+        psLogMsg ("psastro", 3, "failed to find a solution\n");
+        psFree (updates);
+        continue;
+    }
+
+    pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL);
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE | PS_DATA_METADATA, "psastro header stats", updates);
+    psFree (updates);
+  }
+
+  psFree (map);
   psFree (view);
   return true;