Index: /trunk/psastro/Makefile =================================================================== --- /trunk/psastro/Makefile (revision 7013) +++ /trunk/psastro/Makefile (revision 7014) @@ -22,32 +22,24 @@ PSASTRO = \ $(SRC)/psastro.$(ARCH).o \ -$(SRC)/psastroIO.$(ARCH).o \ -$(SRC)/psastroBuildFPA.$(ARCH).o \ $(SRC)/psastroArguments.$(ARCH).o \ -$(SRC)/psastroUtils.$(ARCH).o \ -$(SRC)/psModUtils.$(ARCH).o \ -$(SRC)/psLibUtils.$(ARCH).o \ -$(SRC)/pmAstrom.$(ARCH).o \ -$(SRC)/pmAstromGrid.$(ARCH).o - -PSASTRO-MKTEST = \ -$(SRC)/psastro-mktest.$(ARCH).o \ -$(SRC)/psastroArguments.$(ARCH).o \ -$(SRC)/psastroIO.$(ARCH).o \ -$(SRC)/psModUtils.$(ARCH).o \ -$(SRC)/psLibUtils.$(ARCH).o \ -$(SRC)/pmAstrom.$(ARCH).o +$(SRC)/psastroParseCamera.$(ARCH).o\ +$(SRC)/psastroDataLoad.$(ARCH).o \ +$(SRC)/psastroDataSave.$(ARCH).o \ +$(SRC)/psastroAstromGuess.$(ARCH).o\ +$(SRC)/psastroLoadReferences.$(ARCH).o\ +$(SRC)/psastroConvert.$(ARCH).o \ +$(SRC)/psastroChipAstrom.$(ARCH).o\ +$(SRC)/psastroOneChip.$(ARCH).o\ +$(SRC)/psastroUtils.$(ARCH).o\ +$(SRC)/psastroTestFuncs.$(ARCH).o\ +$(SRC)/psastroWCS.$(ARCH).o psastro: $(BIN)/psastro.$(ARCH) $(BIN)/psastro.$(ARCH) : $(PSASTRO) -$(PSASTRO): $(SRC)/psastro.h $(SRC)/pmAstrom.h +$(PSASTRO): $(SRC)/psastro.h -psastro-mktest: $(BIN)/psastro-mktest.$(ARCH) -$(BIN)/psastro-mktest.$(ARCH) : $(PSASTRO-MKTEST) -$(PSASTRO-MKTEST): $(SRC)/psastro.h $(SRC)/pmAstrom.h +INSTALL = psastro -INSTALL = psastro psastro-mktest - -all: psastro.install psastro-mktest.install +all: psastro.install # dependancy rules for binary code ######################### Index: unk/psastro/src/pmAstrom.c =================================================================== --- /trunk/psastro/src/pmAstrom.c (revision 7013) +++ (revision ) @@ -1,223 +1,0 @@ -# include "pmAstrom.h" - -// sort by mag (descending) -int pmAstromObjSortByFPX (const void **a, const void **b) -{ - pmAstromObj *A = *(pmAstromObj **)a; - pmAstromObj *B = *(pmAstromObj **)b; - - psF32 diff = A->FP.x - B->FP.x; - if (diff > FLT_EPSILON) return (+1); - if (diff < FLT_EPSILON) return (-1); - return (0); -} - -// sort by mag (descending) -int pmAstromObjSortByMag (const void **a, const void **b) -{ - pmAstromObj *A = *(pmAstromObj **)a; - pmAstromObj *B = *(pmAstromObj **)b; - - psF32 diff = A->Mag - B->Mag; - if (diff > FLT_EPSILON) return (-1); - if (diff < FLT_EPSILON) return (+1); - return (0); -} - -psArray *pmAstromRadiusMatch (psArray *st1, psArray *st2, psMetadata *config) { - - bool status; - int jStart; - double dX, dY, dR; - - double RADIUS = psMetadataLookupF32 (&status, config, "PSASTRO.MATCH.RADIUS"); - double RADIUS_SQR = PS_SQR (RADIUS); - - // sort both lists by Focal Plane X coord - st1 = psArraySort (st1, pmAstromObjSortByFPX); - st2 = psArraySort (st2, pmAstromObjSortByFPX); - - psArray *matches = psArrayAlloc (100); - matches->n = 0; - - int i = 0; - int j = 0; - while ((i < st1->n) && (j < st2->n)) { - - dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x; - if (dX < -RADIUS) { - i++; - continue; - } - if (dX > +RADIUS) { - j++; - continue; - } - - jStart = j; - while ((dX > -RADIUS) && (j < st2->n)) { - - dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x; - dY = ((pmAstromObj *)st1->data[i])->FP.y - ((pmAstromObj *)st2->data[j])->FP.y; - dR = dX*dX + dY*dY; - - if (dR > RADIUS_SQR) { - j++; - continue; - } - - // got a match; add to output list - pmAstromMatch *match = pmAstromMatchAlloc (i, j); - psArrayAdd (matches, 100, match); - - j++; - } - j = jStart; - i++; - } - return (matches); -} - -// take two matched star lists and fit a psPlaneTransform between them -// -psPlaneTransform *pmAstromMatchFit (psPlaneTransform *map, psArray *raw, psArray *ref, psArray *match, psMetadata *config) { - - bool status; - pmAstromObj *rawStar, *refStar; - pmAstromMatch *pair; - - if (map == NULL) { - int nX = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NX"); - if (!status) nX = 1; - int nY = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NY"); - if (!status) nY = 1; - map = psPlaneTransformAlloc (nX, nY); - // XXX EAM : not clear that we are allowed to use orders > 1 - } - - psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - stats->clipSigma = psMetadataLookupF32 (&status, config, "PSASTRO.CHIP.NSIGMA"); - stats->clipIter = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NITER"); - - // XXX EAM : clip fit seems to only work for F32! - // XXX EAM : clip fit fails to handle NULL error - psVector *X = psVectorAlloc (match->n, PS_TYPE_F32); - psVector *Y = psVectorAlloc (match->n, PS_TYPE_F32); - psVector *x = psVectorAlloc (match->n, PS_TYPE_F32); - psVector *y = psVectorAlloc (match->n, PS_TYPE_F32); - psVector *wt = psVectorAlloc (match->n, PS_TYPE_F32); - - // take the matched stars, first fit - for (int i = 0; i < match->n; i++) { - - pair = match->data[i]; - rawStar = raw->data[pair->i1]; - refStar = ref->data[pair->i2]; - - X->data.F32[i] = rawStar->chip.x; - Y->data.F32[i] = rawStar->chip.y; - - x->data.F32[i] = refStar->FP.x; - y->data.F32[i] = refStar->FP.y; - - wt->data.F32[i] = 1.0; - } - - // no masking, no errors - psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, wt, x, y); - psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, wt, x, y); - psFree (x); - psFree (y); - psFree (X); - psFree (Y); - psFree (wt); - psFree (stats); - - return (map); -} - -// rotate the focal-plane coordinates about the center coordinate -// angle specified in radians -psArray *pmAstromRotateObj (psArray *old, psPlane center, double angle) { - - double X, Y; - pmAstromObj *newObj; - pmAstromObj *oldObj; - - psArray *new = psArrayAlloc (old->n); - - double cs = cos(angle); - double sn = sin(angle); - double xCenter = center.x; - double yCenter = center.y; - - for (int i = 0; i < old->n; i++) { - - oldObj = (pmAstromObj *)old->data[i]; - newObj = pmAstromObjCopy (oldObj); - - X = oldObj->FP.x - xCenter; - Y = oldObj->FP.y - yCenter; - - newObj->FP.x = X*cs + Y*sn + xCenter; - newObj->FP.y = Y*cs - X*sn + yCenter; - - new->data[i] = newObj; - } - return (new); -} - -static void pmAstromObjFree (pmAstromObj *obj) { - - if (obj == NULL) return; - return; -} - -pmAstromObj *pmAstromObjAlloc (void) { - - pmAstromObj *obj = psAlloc (sizeof(pmAstromObj)); - psMemSetDeallocator (obj, (psFreeFunc) pmAstromObjFree); - - obj->pix.x = 0; - obj->pix.y = 0; - - obj->FP.x = 0; - obj->FP.y = 0; - - obj->TP.x = 0; - obj->TP.y = 0; - - obj->sky.r = 0; - obj->sky.d = 0; - - obj->Mag = 0; - obj->dMag = 0; - - return (obj); -} - -pmAstromObj *pmAstromObjCopy (pmAstromObj *old) { - - pmAstromObj *obj = pmAstromObjAlloc (); - - obj[0] = old[0]; - - return (obj); -} - -static void pmAstromMatchFree (pmAstromMatch *match) { - - if (match == NULL) return; - return; -} - -pmAstromMatch *pmAstromMatchAlloc (int i1, int i2) { - - pmAstromMatch *match = psAlloc (sizeof(pmAstromMatch)); - psMemSetDeallocator(match, (psFreeFunc) pmAstromMatchFree); - - match->i1 = i1; - match->i2 = i2; - - return (match); -} Index: unk/psastro/src/pmAstrom.h =================================================================== --- /trunk/psastro/src/pmAstrom.h (revision 7013) +++ (revision ) @@ -1,44 +1,0 @@ -# include -# include // for strcasecmp -# include // for unlink -# include -# include - -typedef struct { - psPlane pix; - psPlane cell; - psPlane chip; - psPlane FP; - psPlane TP; - psSphere sky; - double Mag, dMag; -} pmAstromObj; - -typedef struct { - int i1; - int i2; -} pmAstromMatch; - -typedef struct { - psPlane center; - psPlane offset; - double angle; - double minMetric; - double minVar; - int nMatch; -} pmAstromStats; - -/* in pmAstrom.c */ -pmAstromObj *pmAstromObjAlloc (void); -pmAstromObj *pmAstromObjCopy (pmAstromObj *old); -pmAstromMatch *pmAstromMatchAlloc (int i1, int i2); -psArray *pmAstromRadiusMatch (psArray *st1, psArray *st2, psMetadata *config); -psArray *pmAstromRotateObj (psArray *old, psPlane center, double angle); -psPlaneTransform *pmAstromMatchFit (psPlaneTransform *map, psArray *st1, psArray *st2, psArray *match, psMetadata *config); -int pmAstromObjSortByFPX (const void **a, const void **b); -int pmAstromObjSortByMag (const void **a, const void **b); - -/* in pmAstromGrid.c */ -pmAstromStats pmAstromGridAngle (psArray *st1, psArray *st2, psMetadata *config); -psPlaneTransform *pmAstromGridApply (psPlaneTransform *map, pmAstromStats stat); -pmAstromStats pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config); Index: unk/psastro/src/pmAstromGrid.c =================================================================== --- /trunk/psastro/src/pmAstromGrid.c (revision 7013) +++ (revision ) @@ -1,218 +1,0 @@ -# include "pmAstrom.h" - -static double maxOffpix; // maximum allowed offset between lists, in raw pixels -static double Scale; // grid pixel scale -static double Offset; // deltas to pixels - -// local function to convert x,y coords to grid bins -// it requires the globals defined above -static bool AstromGridBin (int *dx, int *dy, double dX, double dY) { - - if (fabs(dX) > maxOffpix) return false; - if (fabs(dY) > maxOffpix) return false; - - *dx = dX / Scale + Offset; - *dy = dY / Scale + Offset; - return true; -} - -// match the two lists using the binned delta-delta max -pmAstromStats pmAstromGridAngle (psArray *raw, psArray *ref, psMetadata *config) { - - bool status; - int nPix; // size of matching grid - int nPixHalf; // half-size of matching grid - double dX, dY; // offset between a possible matched pair - int iX, iY; // corresponding grid bin - - pmAstromObj *ob1, *ob2; // short-cut pointers to the objects - pmAstromStats stats; // output match statistics - - // max allowed offset in either X or Y directions - double gridOffset = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.OFFSET"); - - // sampling scale of the grid - double gridScale = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.SCALE"); - - // set the static scaling factors - nPixHalf = (int)(gridOffset / gridScale + 0.5); // half-grid - nPix = 2*nPixHalf + 1; // full grid width - - // these are globals used by p_pmAstromGridBin - maxOffpix = gridScale * (nPixHalf + 0.5); // max offset from true center - Offset = maxOffpix / gridScale; - Scale = gridScale; - - // images used as accumulators for the loop below - psImage *gridNP = psImageAlloc (nPix, nPix, PS_TYPE_U32); - psImage *gridDX = psImageAlloc (nPix, nPix, PS_TYPE_F32); - psImage *gridDY = psImageAlloc (nPix, nPix, PS_TYPE_F32); - psImage *gridD2 = psImageAlloc (nPix, nPix, PS_TYPE_F32); - psImageInit (gridNP, 0); - psImageInit (gridDX, 0); - psImageInit (gridDY, 0); - psImageInit (gridD2, 0); - - // short-cut names for grid images - psS32 **NP = gridNP->data.S32; - psF32 **DX = gridDX->data.F32; - psF32 **DY = gridDY->data.F32; - psF32 **D2 = gridD2->data.F32; - - // accumulate grids for focal plane (L,M) matches - for (int i = 0; i < raw->n; i++) { - ob1 = (pmAstromObj *)raw->data[i]; - for (int j = 0; j < ref->n; j++) { - ob2 = (pmAstromObj *)ref->data[j]; - dX = ob1->FP.x - ob2->FP.x; - dY = ob1->FP.y - ob2->FP.y; - - // fprintf (stderr, "dX,dY: %8.2f %8.2f : %8.2f %8.2f : %8.2f %8.2f\n", dX, dY, ob1->FP.x, ob2->FP.x, ob1->FP.y, ob2->FP.y); - // find bin coordinates for this delta-delta - if (!AstromGridBin (&iX, &iY, dX, dY)) { - continue; // matched pair is too far offset - } - - // accumulate bin stats - NP[iY][iX] ++; - DX[iY][iX] += dX; - DY[iY][iX] += dY; - D2[iY][iX] += PS_SQR(dX) + PS_SQR(dY); - } - } - - // now assess the grid images - { - double minMetric = 1e10; - double minVar = 1e10; - int minX = -1; - int minY = -1; - double metric, var; - - // find the max pixel - psStats *imStats = psStatsAlloc (PS_STAT_MAX); - imStats = psImageStats (imStats, gridNP, NULL, 0); - - // only check bins with at least 1/2 of max bin - int minNpts = 0.5*imStats->max; - fprintf (stderr, "minNpts: %d, max: %d\n", minNpts, (int)(imStats->max)); - - // find the 'best' bin - for (int j = 0; j < gridNP->numRows; j++) { - for (int i = 0; i < gridNP->numCols; i++) { - - if (NP[j][i] < minNpts) continue; - - // this metric emphasize a narrow peak with lots of sources over one with few. - var = fabs((D2[j][i]/NP[j][i]) - PS_SQR(DX[j][i]/NP[j][i]) - PS_SQR(DY[j][i]/NP[j][i])); - metric = var / PS_SQR(PS_SQR(NP[j][i])); - - fprintf (stderr, "try : %f %f (%d pts, %f var, %f met)\n", DX[j][i]/NP[j][i], DY[j][i]/NP[j][i], NP[j][i], var, metric); - - if (metric < minMetric) { - minMetric = metric; - minVar = var; - minX = i; - minY = j; - } - } - } - - // convert the bin to delta-delta - stats.offset.x = DX[minY][minX] / NP[minY][minX]; - stats.offset.y = DY[minY][minX] / NP[minY][minX]; - stats.minMetric = minMetric; - stats.minVar = minVar; - stats.nMatch = NP[minY][minX]; - - // XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better - } - return (stats); -} - -// match two star lists -pmAstromStats pmAstromGridMatch (psArray *raw, psArray *ref, psMetadata *config) { - - bool status; - double xMin, xMax, yMin, yMax; - pmAstromObj *obj; - psArray *rot; - - pmAstromStats minStat, newStat; - psPlane center; - - // find center of the raw field (focal-plane coords) - xMin = yMin = +1e10; - xMax = yMax = -1e10; - for (int i = 0; i < raw->n; i++) { - obj = (pmAstromObj *)raw->data[i]; - xMin = PS_MIN (obj->FP.x, xMin); - xMax = PS_MAX (obj->FP.x, xMax); - yMin = PS_MIN (obj->FP.y, yMin); - yMax = PS_MAX (obj->FP.y, yMax); - } - center.x = 0.5*(xMin + xMax); - center.y = 0.5*(yMin + yMax); - - double minAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.ANGLE"); - double maxAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.ANGLE"); - double delAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.DEL.ANGLE"); - - minStat.minMetric = 1e10; - for (double angle = minAngle; angle <= maxAngle; angle += delAngle) { - rot = pmAstromRotateObj (raw, center, angle); - newStat = pmAstromGridAngle (rot, ref, config); - newStat.angle = angle; - newStat.center = center; - if (newStat.minMetric < minStat.minMetric) { - minStat = newStat; - } - psFree (rot); - } - fprintf (stderr, "best: %f %f (%d pts, %f var, %f met)\n", newStat.offset.x, newStat.offset.y, newStat.nMatch, newStat.minVar, newStat.minMetric); - return (minStat); -} - -// apply the measured FPA offset and rotation (stat) to the fpa astrom structures -psPlaneTransform *pmAstromGridApply (psPlaneTransform *map, pmAstromStats stat) { - - double cs = cos (stat.angle); - double sn = sin (stat.angle); - - double dx = (map->x->coeff[0][0] - stat.center.x); - double dy = (map->y->coeff[0][0] - stat.center.y); - - // new offset - map->x->coeff[0][0] = cs*dx + sn*dy - stat.offset.x + stat.center.x; - map->y->coeff[0][0] = -sn*dx + cs*dy - stat.offset.y + stat.center.y; - - // original rotation matrix - double pc1_1 = map->x->coeff[1][0]; - double pc1_2 = map->x->coeff[0][1]; - double pc2_1 = map->y->coeff[1][0]; - double pc2_2 = map->y->coeff[0][1]; - - // new rotation matrix - map->x->coeff[1][0] = +cs*pc1_1 + sn*pc2_1; - map->x->coeff[0][1] = +cs*pc1_2 + sn*pc2_2; - map->y->coeff[1][0] = -sn*pc1_1 + cs*pc2_1; - map->y->coeff[0][1] = -sn*pc1_2 + cs*pc2_2; - - return (map); -} - -/* Illustration of the grid bins - dX bin - -35:-25 -> 0 bin = dX / Scale + Offset - -25:-15 -> 1 Scale = 10 - -15:-05 -> 2 Offset = 3.5 - -05:+05 -> 3 nPix = 3 (maxOffset = 35 = (nPix + 0.5)*dXix - +05:+15 -> 4 dPix = 10 - +15:+25 -> 5 - +25:+35 -> 6 - - maxOffsetRequest = 30 - nPix = (int) (maxOffset / dPix + 0.5); - maxOffset = (nPix + 0.5)*Scale; -*/ - Index: unk/psastro/src/pmFPAReadSMP.c =================================================================== --- /trunk/psastro/src/pmFPAReadSMP.c (revision 7013) +++ (revision ) @@ -1,134 +1,0 @@ - -typedef struct { - float x; - float y; - float fitMag; - float apMag; - float galMag; - float dMag; - float sky; - float Sx; - float Sy; - float Theta; -} pmStar; - -// Given a FITS file pointer, read the table of object data -bool pmFPAviewReadSMP (pmFPAview *view, pmFPAfile *file) { - - pmFPA *fpa = view->fpa; - - if (view->chip == -1) { - pmFPARead (fpa, view, file); - return true; - } - - if (view->chip >= fpa->chips->n) { - return false; - } - pmChip *chip = fpa->chips->data[view->chip]; - - if (view->cell == -1) { - pmChipRead (chip, view, file); - return true; - } - - if (view->cell >= chip->cells->n) { - return false; - } - pmCell *cell = chip->cells->data[view->cell]; - - if (view->readout == -1) { - pmCellRead (cell, view, file); - return true; - } - - if (view->readout >= cell->readouts->n) { - return false; - } - pmReadout *readout = cell->readouts->data[view->readout]; - - pmReadoutReadSMP (readout, view, file); - return true; -} - -// read in all chip-level SMP files for this FPA -pmFPAReadSMP (pmFPA *fpa, pmFPAview *view, pmFPAfile *file) { - - for (int i = 0; i < fpa->chips->n; i++) { - - pmChip *chip = fpa->data[i]; - pmChipReadSMP (chip, view, file); - } - return true; -} - -// read in all cell-level SMP files for this chip -pmChipReadSMP (pmChip *chip, pmFPAview *view, pmFPAfile *file) { - - for (int i = 0; i < chip->cells->n; i++) { - - pmCell *cell = chip->data[i]; - pmCellReadSMP (cell, view, file); - } - return true; -} - -// read in all readout-level SMP files for this cell -pmCellReadSMP (pmCell *cell, pmFPAview *view, pmFPAfile *file) { - - for (int i = 0; i < cell->readouts->n; i++) { - - pmReadout *readout = cell->data[i]; - pmReadoutReadSMP (readout, view, file); - } - return true; -} - -// read in a readout from the fits file -pmReadoutReadSMP (pmReadout *readout, pmFPAview *view, pmFPAfile *file) { - - psFits *fits = file->fits; - - // there is one SMP header per cell and one SMP table per readout - psMetadata *header = pmReadoutGetHeader (readout); - if (header == NULL) { - // determine the extname for the HEAD - headname = pmConfigNameFromRule (file->extextra, camera, view); - psFitsMoveExtName(fits, headname); - psMetadata *header = psFitsReadHeader (NULL, fits); - } - - // determine the extname for the HEAD - dataname = pmConfigNameFromRule (file->extname, camera, view); - - // find config information for output header - float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); - - psFitsMoveExtName(fits, dataname); - psArray *table = psFitsReadTable (fits); - - // validate a single row of the table (must match SMP) - - psArray *sources = psArrayAlloc (table->n); - - // convert the table to the pmSource entries - for (int i = 0; i < table->n; i++) { - pmStar *source = pmSourceAlloc (); - - psMetadata *row = table->data[i]; - - source->x = psMetadataLookupF32 (&status, row, "X_PIX"); - source->y = psMetadataLookupF32 (&status, row, "Y_PIX"); - source->fitMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT; - source->galMag = psMetadataLookupF32 (&status, row, "MAG_GAL") - ZERO_POINT; - source->apMag = psMetadataLookupF32 (&status, row, "MAG_AP") - ZERO_POINT; - source->dMag = psMetadataLookupF32 (&status, row, "MAG_ERR"); - source->sky = psMetadataLookupF32 (&status, row, "LOG_SKY"); - source->sx = psMetadataLookupF32 (&status, row, "FWHM_X"); - source->sy = psMetadataLookupF32 (&status, row, "FWHM_Y"); - source->theta = psMetadataLookupF32 (&status, row, "THETA"); - - sources->data[i] = source; - } - status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources); -} Index: unk/psastro/src/psLibUtils.c =================================================================== --- /trunk/psastro/src/psLibUtils.c (revision 7013) +++ (revision ) @@ -1,204 +1,0 @@ -# include // for strncasecmp -# include -# include "psLibUtils.h" - -// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae) -// 2 = default, -v = 3, -vv = 4, -vvv = 5 -psS32 psLogArguments (int *argc, char **argv) { - - int N, level; - - // default log level is 2 - level = 2; - - // set in order, so that -vvv overrides -vv overrides -v - if ((N = psArgumentGet (*argc, argv, "-v"))) { - psArgumentRemove (N, argc, argv); - level = 3; - } - if ((N = psArgumentGet (*argc, argv, "-vv"))) { - psArgumentRemove (N, argc, argv); - level = 4; - } - if ((N = psArgumentGet (*argc, argv, "-vvv"))) { - psArgumentRemove (N, argc, argv); - level = 5; - } - - if ((N = psArgumentGet (*argc, argv, "-logfmt"))) { - if (*argc < N + 2) { - psAbort ("psLogArguments", "USAGE: -logfmt (format)"); - } - psArgumentRemove (N, argc, argv); - psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid - psArgumentRemove (N, argc, argv); - } - - // set the level, return the level - psLogSetLevel (level); - return (level); -} - -// set trace levels by facility -psS32 psTraceArguments (int *argc, char **argv) { - - int N; - - // argument format is: -trace (facil) (level) - while ((N = psArgumentGet (*argc, argv, "-trace"))) { - if (*argc < N + 3) { - psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)"); - } - psArgumentRemove (N, argc, argv); - psTraceSetLevel (argv[N], atoi(argv[N+1])); - psArgumentRemove (N, argc, argv); - psArgumentRemove (N, argc, argv); - } - if ((N = psArgumentGet (*argc, argv, "-trace-levels"))) { - psTracePrintLevels (); - exit (2); - } - return (TRUE); -} - -psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input) { - - psPolynomial2D *output = psPolynomial2DAlloc (input->nX, input->nY, input->type); - - for (int i = 0; i < input->nX; i++) { - for (int j = 0; j < input->nY; j++) { - output->mask[i][j] = input->mask[i][j]; - output->coeff[i][j] = input->coeff[i][j]; - output->coeffErr[i][j] = input->coeffErr[i][j]; - } - } - return (output); -} - -psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input) { - - psPolynomial4D *output = psPolynomial4DAlloc (input->nX, input->nY, input->nZ, input->nT, input->type); - - for (int i = 0; i < input->nX; i++) { - for (int j = 0; j < input->nY; j++) { - for (int k = 0; k < input->nZ; k++) { - for (int m = 0; m < input->nT; m++) { - output->mask[i][j][k][m] = input->mask[i][j][k][m]; - output->coeff[i][j][k][m] = input->coeff[i][j][k][m]; - output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m]; - } - } - } - } - return (output); -} - -psPlaneDistort *psPlaneDistortCopy (psPlaneDistort *input) { - - psPlaneDistort *output = psPlaneDistortAlloc (input->x->nX, input->x->nY, input->x->nZ, input->x->nT); - - for (int i = 0; i < input->x->nX; i++) { - for (int j = 0; j < input->x->nY; j++) { - for (int k = 0; k < input->x->nZ; k++) { - for (int m = 0; m < input->x->nT; m++) { - // x-terms - output->x->mask[i][j][k][m] = input->x->mask[i][j][k][m]; - output->x->coeff[i][j][k][m] = input->x->coeff[i][j][k][m]; - output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m]; - // y-terms - output->y->mask[i][j][k][m] = input->y->mask[i][j][k][m]; - output->y->coeff[i][j][k][m] = input->y->coeff[i][j][k][m]; - output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m]; - } - } - } - } - return (output); -} - -psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input) { - - psPlaneTransform *output = psPlaneTransformAlloc (input->x->nX, input->x->nY); - - for (int i = 0; i < input->x->nX; i++) { - for (int j = 0; j < input->x->nY; j++) { - // x-terms - output->x->mask[i][j] = input->x->mask[i][j]; - output->x->coeff[i][j] = input->x->coeff[i][j]; - output->x->coeffErr[i][j] = input->x->coeffErr[i][j]; - // y-terms - output->y->mask[i][j] = input->y->mask[i][j]; - output->y->coeff[i][j] = input->y->coeff[i][j]; - output->y->coeffErr[i][j] = input->y->coeffErr[i][j]; - } - } - return (output); -} - -psProjection *psProjectionCopy (psProjection *input) { - - psProjection *output = psProjectionAlloc (input->R, input->D, input->Xs, input->Ys, input->type); - return (output); -} - -// very crude distortion inversion: assumes 0 order in z and t, linear in x and y: -psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort) { - PS_ASSERT_PTR_NON_NULL(distort, 0); - PS_ASSERT_PTR_NON_NULL(distort->x, 0); - PS_ASSERT_PTR_NON_NULL(distort->y, 0); - - psPlaneDistort *out = psPlaneDistortAlloc(1, 1, 0, 0); - - /* simple matrix inversion code */ - - psF64 r11 = distort->x->coeff[1][0][0][0]; - psF64 r12 = distort->x->coeff[0][1][0][0]; - psF64 r21 = distort->y->coeff[1][0][0][0]; - psF64 r22 = distort->y->coeff[0][1][0][0]; - psF64 xo = distort->x->coeff[0][0][0][0]; - psF64 yo = distort->y->coeff[0][0][0][0]; - - psF64 invDet = 1.0 / (r11 * r22 - r12 * r21); // Inverse of the determinant - - out->x->coeff[1][0][0][0] = +invDet * r22; - out->x->coeff[0][1][0][0] = -invDet * r12; - out->y->coeff[1][0][0][0] = -invDet * r21; - out->y->coeff[0][1][0][0] = +invDet * r11; - - out->x->coeff[0][0][0][0] = - invDet * (r22 * xo - r12 * yo); - out->y->coeff[0][0][0][0] = - invDet * (r11 * yo - r21 * xo); - - return(out); -} - -// returns the rotation term, forcing positive parity -double psPlaneTransformGetRotation (psPlaneTransform *map) { - - if (map->x->nX < 1) return 0; - if (map->x->nY < 1) return 0; - - if (map->y->nX < 1) return 0; - if (map->y->nY < 1) return 0; - - double pc1_1 = map->x->coeff[1][0]; - double pc1_2 = map->x->coeff[0][1]; - double pc2_1 = map->y->coeff[1][0]; - double pc2_2 = map->y->coeff[0][1]; - - double px = SIGN (pc1_1); - double py = SIGN (pc2_2); - - // both x and y terms imply an angle. take the average - double t1 = -atan2 (px*pc1_2, px*pc1_1); - double t2 = +atan2 (py*pc2_1, py*pc2_2); - - // careful near -pi,+pi boundary... - if (t1 - t2 > M_PI/2) t2 += 2*M_PI; - if (t2 - t1 > M_PI/2) t1 += 2*M_PI; - - double theta = 0.5*(t1 + t2); - while (theta < M_PI) theta += 2*M_PI; - while (theta > M_PI) theta -= 2*M_PI; - - return (theta); -} Index: unk/psastro/src/psLibUtils.h =================================================================== --- /trunk/psastro/src/psLibUtils.h (revision 7013) +++ (revision ) @@ -1,22 +1,0 @@ -# ifndef PS_LIB_UTILS -# define PS_LIB_UTILS - -# define psMemCopy(A)(psMemIncrRefCounter((A))) -# define DEG_RAD 57.295779513082322 -# define RAD_DEG 0.017453292519943 - -# define SIGN(X) (((X) == 0) ? 0 : ((fabs((double)(X))) / (X))) - -// psLib extra utilities -psS32 psLogArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity) -psS32 psTraceArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity) - -psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input); -psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input); -psPlaneDistort *psPlaneDistortCopy (psPlaneDistort *input); -psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input); -psProjection *psProjectionCopy (psProjection *input); -psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort); -double psPlaneTransformGetRotation (psPlaneTransform *map); - -# endif Index: unk/psastro/src/psastro-mktest.c =================================================================== --- /trunk/psastro/src/psastro-mktest.c (revision 7013) +++ (revision ) @@ -1,169 +1,0 @@ -# include "psastro.h" - -int main (int argc, char **argv) { - - bool status; - - // load configuration information - psMetadata *config = testArguments (&argc, argv); - - // fill in the header with astrometry info as we go... - psMetadata *header = psMetadataAlloc(); - - // make fake star set - psRandom *rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); - - int NX = psMetadataLookupS32 (&status, config, "NX"); - int NY = psMetadataLookupS32 (&status, config, "NY"); - int NSTARS = psMetadataLookupS32 (&status, config, "NSTARS"); - - psMetadataAdd (header, PS_LIST_TAIL, "SIMPLE", PS_DATA_BOOL | PS_META_REPLACE, "foo", true); - psMetadataAdd (header, PS_LIST_TAIL, "BITPIX", PS_DATA_S32 | PS_META_REPLACE, "foo", 16); - psMetadataAdd (header, PS_LIST_TAIL, "NAXIS", PS_DATA_S32 | PS_META_REPLACE, "foo", 2); - psMetadataAdd (header, PS_LIST_TAIL, "NAXIS1", PS_DATA_S32 | PS_META_REPLACE, "foo", NX); - psMetadataAdd (header, PS_LIST_TAIL, "NAXIS2", PS_DATA_S32 | PS_META_REPLACE, "foo", NY); - - psArray *rawstars = psArrayAlloc (NSTARS); - for (int i = 0; i < NSTARS; i++) { - pmAstromObj *star = pmAstromObjAlloc (); - star->pix.x = NX * psRandomUniform (rnd); - star->pix.y = NY * psRandomUniform (rnd); - star->Mag = 10.0; - star->dMag = 0.05; - rawstars->data[i] = star; - } - - pmFPA *fpa = pmFPAAlloc (NULL, NULL); - pmChip *chip = pmChipAlloc (fpa); - pmCell *cell = pmCellAlloc (chip); - pmReadout *readout = pmReadoutAlloc (cell); - - { // define readout -> cell transformation - readout->col0 = psMetadataLookupF32 (&status, config, "READOUT.COL.0"); - readout->row0 = psMetadataLookupF32 (&status, config, "READOUT.ROW.0"); - readout->colBins = psMetadataLookupF32 (&status, config, "READOUT.COL.BIN"); - readout->rowBins = psMetadataLookupF32 (&status, config, "READOUT.ROW.BIN"); - } - - { // define cell -> chip transformation - cell->toChip = psPlaneTransformAlloc (1, 1); - cell->toChip->x->coeff[0][0] = psMetadataLookupF32 (&status, config, "CELL.COL.0"); - cell->toChip->y->coeff[0][0] = psMetadataLookupF32 (&status, config, "CELL.ROW.0"); - - cell->toChip->x->coeff[1][0] = 1; - cell->toChip->x->coeff[0][1] = 0; - cell->toChip->x->coeff[1][1] = 0; - - cell->toChip->y->coeff[1][0] = 0; - cell->toChip->y->coeff[0][1] = 1; - cell->toChip->y->coeff[1][1] = 0; - } - - { // define chip -> fpa transformation - chip->toFPA = psPlaneTransformAlloc (1, 1); - chip->toFPA->x->coeff[0][0] = psMetadataLookupF32 (&status, config, "CHIP.COL.0"); - chip->toFPA->y->coeff[0][0] = psMetadataLookupF32 (&status, config, "CHIP.ROW.0"); - - double theta = psMetadataLookupF32 (&status, config, "CHIP.THETA"); - double px = psMetadataLookupF32 (&status, config, "CHIP.PARITY.X"); - double py = psMetadataLookupF32 (&status, config, "CHIP.PARITY.Y"); - - double cs = cos (theta*RAD_DEG); - double sn = sin (theta*RAD_DEG); - - chip->toFPA->x->coeff[1][0] = +px*cs; - chip->toFPA->x->coeff[0][1] = +px*sn; - chip->toFPA->x->coeff[1][1] = 0; - - chip->toFPA->y->coeff[1][0] = -py*sn; - chip->toFPA->y->coeff[0][1] = +py*cs; - chip->toFPA->y->coeff[1][1] = 0; - - chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); - } - - { // define fpa -> tpa transformation - fpa->toTangentPlane = psPlaneDistortAlloc (1, 1, 0, 0); - fpa->toTangentPlane->x->coeff[0][0][0][0] = psMetadataLookupF32 (&status, config, "FPA.COL.0"); - fpa->toTangentPlane->y->coeff[0][0][0][0] = psMetadataLookupF32 (&status, config, "FPA.ROW.0"); - - double theta = psMetadataLookupF32 (&status, config, "FPA.THETA"); - // double px = psMetadataLookupF32 (&status, config, "FPA.PARITY.X"); - // double py = psMetadataLookupF32 (&status, config, "FPA.PARITY.Y"); - - double cs = cos (theta*RAD_DEG); - double sn = sin (theta*RAD_DEG); - - fpa->toTangentPlane->x->coeff[1][0][0][0] = +cs; - fpa->toTangentPlane->x->coeff[0][1][0][0] = +sn; - fpa->toTangentPlane->x->coeff[1][1][0][0] = 0; - - fpa->toTangentPlane->y->coeff[1][0][0][0] = -sn; - fpa->toTangentPlane->y->coeff[0][1][0][0] = +cs; - fpa->toTangentPlane->y->coeff[1][1][0][0] = 0; - - fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane); - } - - { // define tpa -> sky projection, update output header - - double RA = psMetadataLookupF32 (&status, config, "RA"); - double DEC = psMetadataLookupF32 (&status, config, "DEC"); - double PS = psMetadataLookupF32 (&status, config, "PLATE.SCALE"); - // XXX EAM : PS units are inverted in psDeproject (fixed in eam branch) - PS /= 3600.0 * DEG_RAD; // convert from arcsec to radians - RA *= RAD_DEG; // convert from degrees to radians - DEC *= RAD_DEG; // convert from degrees to radians - - fpa->projection = psProjectionAlloc (RA, DEC, PS, PS, PS_PROJ_SIN); - - // find the readout coord for ref pixel: - psSphere tmpSky; - psPlane tmp1, tmp2; - tmpSky.r = RA; - tmpSky.d = DEC; - - psCoordSkyToTP (&tmp1, &tmpSky, fpa->projection); - psCoordTPToFPA (&tmp2, &tmp1, 0.0, 0.0, fpa); - psCoordFPAToChip (&tmp1, &tmp2, chip); - psCoordChipToCell_EAM (&tmp2, &tmp1, cell); - psCoordCellToReadout (&tmp1, &tmp2, readout); - - // write output header WCS info - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "foo", "RA---SIN"); - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "foo", "DEC--SIN"); - psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32 | PS_META_REPLACE, "foo", RA*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32 | PS_META_REPLACE, "foo", DEC*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32 | PS_META_REPLACE, "foo", tmp1.x); - psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32 | PS_META_REPLACE, "foo", tmp1.y); - psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32 | PS_META_REPLACE, "foo", PS*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32 | PS_META_REPLACE, "foo", PS*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[1][0]); - psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[0][1]); - psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[1][0]); - psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[0][1]); - } - - // project stars from readout up to sky - psastroProjectRawstars (rawstars, readout); - - // generate refstars from rawstars (for test comparison) - psArray *refstars = psArrayAlloc (rawstars->n); - for (int i = 0; i < rawstars->n; i++) { - pmAstromObj *oStar = pmAstromObjAlloc (); - pmAstromObj *iStar = rawstars->data[i]; - oStar->sky = iStar->sky; - oStar->Mag = iStar->Mag; - oStar->dMag = iStar->dMag; - refstars->data[i] = oStar; - } - - // project stars from sky down to readout - psastroProjectRefstars (refstars, readout); - - testWriteCMP (header, argv[2], rawstars); - testWriteRef (argv[3], rawstars); - testWriteRaw (argv[4], rawstars); - testWriteRaw (argv[5], refstars); - exit (0); -} Index: /trunk/psastro/src/psastro.c =================================================================== --- /trunk/psastro/src/psastro.c (revision 7013) +++ /trunk/psastro/src/psastro.c (revision 7014) @@ -5,44 +5,60 @@ psTimerStart ("complete"); - psMetadata *header = NULL; + // model inits are needed in pmSourceIO + // XXX do we want to do the local model inits a la psphot? + pmModelGroupInit (); // load configuration information pmConfig *config = psastroArguments (&argc, argv); - // load input data (config and images (signal, noise, mask) + // load identify the data sources psastroParseCamera (config); - // perform the astrometric solution - psastroDataLoop (input, config); + // load the raw pixel data (from PSPHOT.SOURCES) + // select subset of stars for astrometry + psastroDataLoad (config); - // select astrometry stars from the psphot sources - // limit fit to bright stars only - psastroSelectAstromStars (config); + // interpret the available initial astrometric information + // apply the initial guess + psastroAstromGuess (config); // load the reference stars overlapping the data stars - psastroLoadReference (config); + psArray *refs = psastroLoadReferences (config); - // fpa and subset point to the same astrometry terms + psastroChipAstrom (config, refs); + + # if 0 + // perform the desired astrometry analysis switch (mode) { case stack: - psastroStackAstrom (config); + psastroStackAstrom (config, refs); + break + case chip: + psastroChipAstrom (config, refs); break case mosaic: - psastroMosaicAstrom (config); - break - case mosaic: - psastroMosaicAstrom (config); + psastroMosaicAstrom (config, refs); break; default: break; } + # endif - // perform the astrometric solution + // write out the results psastroDataSave (config); - psLogMsg ("psphot", 3, "complete psastro run: %f sec\n", psTimerMark ("complete")); + psLogMsg ("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete")); + psFree (refs); psFree (config); - psFree (fpa); + + psTimerStop (); + psMemCheckCorruption (true); + pmModelGroupCleanup (); + psTimeFinalize (); + pmConceptsDone (); + pmConfigDone (); + fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro"); + // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro"); exit (0); Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 7013) +++ /trunk/psastro/src/psastro.h (revision 7014) @@ -3,7 +3,10 @@ # include // for unlink # include -# include -# include "psLibUtils.h" -# include "pmAstrom.h" +# include + +# define psMemCopy(A)(psMemIncrRefCounter((A))) +# define DEG_RAD 57.295779513082322 +# define RAD_DEG 0.017453292519943 +# define SIGN(X) (((X) == 0) ? 0 : ((fabs((double)(X))) / (X))) # define fromTPA fromTangentPlane @@ -11,56 +14,26 @@ # define toSky projection -// How much of the FPA to load at a time -typedef enum { - PP_LOAD_NONE, // Don't load anything - PP_LOAD_FPA, // Load the entire FPA at once - PP_LOAD_CHIP, // Load by chip - PP_LOAD_CELL, // Load by cell -} ppImageLoadDepth; +pmConfig *psastroArguments (int *argc, char **argv); +bool psastroParseCamera (pmConfig *config); +bool psastroDataLoad (pmConfig *config); +bool psastroDataSave (pmConfig *config); -// Configuration data -typedef struct { - psMetadata *site; // The site configuration - psMetadata *camera; // The camera configuration - psMetadata *recipe; // The recipe (i.e., specific setups) - psMetadata *arguments; // Command-line arguments - psMetadata *options; // Command-line recipe options - psDB *database; // Database handle -} ppConfig; +bool psastroConvert (pmReadout *readout, psMetadata *recipe); +psArray *pmSourceToAstromObj (psArray *sources); +bool psastroAstromGuess (pmConfig *config); -// A file to process -typedef struct { - char *filename; // File name - psFits *fits; // The FITS file handle - psMetadata *phu; // The FITS header - pmFPA *fpa; // The FPA, with pixels and extensions -} ppFile; +bool pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale); +bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale); +psPlaneDistort *psPlaneDistortIdentity (); -bool pmAstromReadWCS (psPlaneTransform **toFPA, psProjection **toSky, psMetadata *header); -bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header); +psArray *psastroLoadReferences (pmConfig *config); +bool psastroChipAstrom (pmConfig *config, psArray *refs); +bool psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refset, psArray *rawset, psMetadata *recipe); -pmFPA *pmFPACopyAstrom (pmFPA *inFPA); -psMetadata *psastroArguments (int *argc, char **argv); -pmFPA *psastroBuildFPA (psMetadata *header, psArray *stars); -bool psastroChipAstrom (pmFPA *fpa, psMetadata *config); -psArray *psastroLoadReference (psMetadata *header, psMetadata *config); -psArray *psastroReadCMP (psMetadata **header, char *filename); -bool psastroSelectBrightStars (pmFPA *fpa, psMetadata *config); -bool psastroWriteCMP (pmFPA *fpa, char *filename); +// bool psVectorSmooth (psVector *vector, double sigma, double Nsigma); +// utility functions: +bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars); +bool psastroWriteStars (char *filename, psArray *sources); +bool psastroWriteTransform (psPlaneTransform *map); -psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout); -psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout); -psPlane *psCoordChipToCell_EAM(psPlane* cellCoord, const psPlane* chipCoord, const pmCell* cell); - -bool testWriteCMP (psMetadata *header, char *filename, psArray *sources); -bool testWriteRef (char *filename, psArray *sources); -bool testWriteRaw (char *filename, psArray *sources); -psMetadata *testArguments (int *argc, char **argv); - -bool psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky); -bool psastroProjectRawstars (psArray *stars, pmReadout *readout); -bool psastroProjectRefstars (psArray *stars, pmReadout *readout); - -psPlaneTransform *p_psPlaneTransformLinearInvert_EAM(psPlaneTransform *transform); -pmFPA *pmFPACopyAstrom (pmFPA *inFPA); Index: /trunk/psastro/src/psastroArguments.c =================================================================== --- /trunk/psastro/src/psastroArguments.c (revision 7013) +++ /trunk/psastro/src/psastroArguments.c (revision 7014) @@ -1,4 +1,5 @@ # include "psastro.h" # include +// XXX leak free 2006.04.27 static void usage (void) { @@ -9,4 +10,5 @@ pmConfig *psastroArguments (int *argc, char **argv) { + bool status; int N; @@ -16,5 +18,5 @@ psLogSetFormat ("M"); - // these other options override the PSPHOT recipe options + // these other options override the PSASTRO recipe options psMetadata *options = psMetadataAlloc (); @@ -27,5 +29,5 @@ // load config data from default locations - pmConfig *config = pmConfigRead(argc, argv); + pmConfig *config = pmConfigRead (argc, argv); // Storage for other command-line arguments @@ -43,5 +45,5 @@ } - status = pmConfigFileSetsMD (config->arguments, "INPUT", "-file", "-list"); + status = pmConfigFileSetsMD (config->arguments, argc, argv, "INPUT", "-file", "-list"); if (!status) { usage ();} Index: /trunk/psastro/src/psastroAstromGuess.c =================================================================== --- /trunk/psastro/src/psastroAstromGuess.c (revision 7014) +++ /trunk/psastro/src/psastroAstromGuess.c (revision 7014) @@ -0,0 +1,110 @@ +# include "psastro.h" + +// XXX this function assumes the initial astrometry arrives in the form of +// XXX WCS keywords in the headers corresponding to the chips. +bool psastroAstromGuess (pmConfig *config) { + + bool newFPA = true; + bool status = false; + double RAmin, RAmax, RAminSky, RAmaxSky; + double DECmin, DECmax; + + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO"); + if (!recipe) { + psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n"); + exit(EXIT_FAILURE); + } + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!input) { + psErrorStackPrint(stderr, "Can't find input data!\n"); + exit(EXIT_FAILURE); + } + + double plateScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLATE.SCALE"); + if (!status) plateScale = 1.0; + + pmFPAview *view = pmFPAviewAlloc (0); + pmFPA *fpa = input->fpa; + + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace (__func__, 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + + // read WCS data from the corresponding header + pmHDU *hdu = pmFPAviewThisHDU (view, fpa); + + pmAstromReadWCS (&fpa->projection, &fpa->toTangentPlane, &chip->toFPA, hdu->header, plateScale); + if (newFPA) { + newFPA = false; + RAminSky = fpa->projection->R - M_PI; + RAmaxSky = fpa->projection->R + M_PI; + RAmin = RAmax = fpa->projection->R; + DECmin = DECmax = fpa->projection->D; + } + + // build projection inversions + // XXX region should be set from image data + psRegion region = {0, 0, 2000, 4500}; + fpa->fromTangentPlane = psPlaneDistortIdentity (); + chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, region, 20); + + // apply the new WCS guess data to all of the data in the readouts + // XXX should this go into a different function? this would separate WCS interpretation from application + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace (__func__, 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + + psArray *objects = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS"); + if (objects == NULL) { continue; } + + for (int i = 0; i < objects->n; i++) { + pmAstromObj *object = objects->data[i]; + + psPlaneTransformApply (object->FP, chip->toFPA, object->chip); + psPlaneDistortApply (object->TP, fpa->toTangentPlane, object->FP, 0.0, 0.0); + p_psDeproject (object->sky, object->TP, fpa->projection); + + // rationalize ra to sky range centered on boresite + while (object->sky->r < RAminSky) object->sky->r += M_PI; + while (object->sky->r > RAmaxSky) object->sky->r -= M_PI; + + RAmin = PS_MIN (object->sky->r, RAmin); + RAmax = PS_MAX (object->sky->r, RAmax); + + DECmin = PS_MIN (object->sky->d, DECmin); + DECmax = PS_MAX (object->sky->d, DECmax); + } + } + } + } + + psLogMsg ("psastro", 2, "loaded data from %f,%f - %f,%f\n", RAmin, DECmin, RAmax, DECmax); + + psMetadataAdd (recipe, PS_LIST_TAIL, "RA_MIN", PS_DATA_F32 | PS_META_REPLACE, "", RAmin); + psMetadataAdd (recipe, PS_LIST_TAIL, "RA_MAX", PS_DATA_F32 | PS_META_REPLACE, "", RAmax); + psMetadataAdd (recipe, PS_LIST_TAIL, "DEC_MIN", PS_DATA_F32 | PS_META_REPLACE, "", DECmin); + psMetadataAdd (recipe, PS_LIST_TAIL, "DEC_MAX", PS_DATA_F32 | PS_META_REPLACE, "", DECmax); + + psFree (view); + return true; +} + +/* coordinate frame hierachy + pixels (on a given readout) + cell + chip + FP (focal plane) + TP (tangent plane) + sky (ra, dec) +*/ Index: unk/psastro/src/psastroBuildFPA.c =================================================================== --- /trunk/psastro/src/psastroBuildFPA.c (revision 7013) +++ (revision ) @@ -1,69 +1,0 @@ -# include "psastro.h" - -pmFPA *psastroBuildFPA (psMetadata *header, psArray *stars) { - - // this function constructs a basic template pmFPA structure based on the information in the header - // metadata. it converts the astrometric information in the header (WCS, etc) to the corresponding - // pmFPA/pmChip, etc distortion and projection elements. - - // this implementation is basic, assuming a single readout,chip,cell - // I will not handle the datasec/biassec information to define the readout/cell - // these things are already handled by Paul's code. this function simple looks for the WCS information - // and provides those conversions. - - int Nchips = 1; - int Ncells = 1; - int Nreadouts = 1; - - // allocate the structures - - pmFPA *fpa = pmFPAAlloc (NULL, NULL); - - // identify matrix for fpa to tpa - fpa->toTPA = psPlaneDistortAlloc (1, 1, 0, 0); - fpa->toTPA->x->coeff[1][0][0][0] = 1; - fpa->toTPA->x->mask[1][1][0][0] = 1; - - fpa->toTPA->y->coeff[0][1][0][0] = 1; - fpa->toTPA->y->mask[1][1][0][0] = 1; - fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane); - - psArrayRealloc (fpa->chips, Nchips); - for (int i = 0; i < Nchips; i++) { - pmChip *chip = pmChipAlloc (fpa); - - psArrayRealloc (chip->cells, Ncells); - for (int j = 0; j < Ncells; j++) { - pmCell *cell = pmCellAlloc (chip); - - // XXX EAM : extend this function to take more than one header - cell->header = header; - pmAstromReadWCS (&chip->toFPA, &fpa->toSky, header); - chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); - - // allocate identity matrix for cell to chip - cell->toChip = psPlaneTransformAlloc (1, 1); - cell->toChip->x->coeff[1][0] = 1; - cell->toChip->x->mask[1][1] = 1; - cell->toChip->y->coeff[0][1] = 1; - cell->toChip->y->mask[1][1] = 1; - - psArrayRealloc (cell->readouts, Nreadouts); - for (int k = 0; k < Nreadouts; k++) { - pmReadout *readout = pmReadoutAlloc (cell); - - // XXX EAM : extend this function to take more than one stars... - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.FULLSET", PS_DATA_ARRAY, "stars from analysis", stars); - - // XXX EAM : this information should be in the header... - readout->col0 = readout->row0 = 0; - readout->colBins = readout->rowBins = 1; - - cell->readouts->data[j] = readout; - } - chip->cells->data[j] = cell; - } - fpa->chips->data[i] = chip; - } - return (fpa); -} Index: /trunk/psastro/src/psastroChipAstrom.c =================================================================== --- /trunk/psastro/src/psastroChipAstrom.c (revision 7013) +++ /trunk/psastro/src/psastroChipAstrom.c (revision 7014) @@ -1,73 +1,87 @@ # include "psastro.h" -// measure per-chip astrometry terms -// this function requires a chip from a well-formed pmFPA -// with PSASTRO.OBJECTS on the readouts -bool psastroChipAstrom (pmChip *chip, psMetadata *config) { +bool psastroChipAstrom (pmConfig *config, psArray *refs) { bool status; - psArray *match; - pmAstromStats stats; + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; - // cells->n > 1 is not well-defined - if (chip->cells->n > 1) { - psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nCells > 1"); - return false; + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO"); + if (!recipe) { + psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n"); + exit(EXIT_FAILURE); } - pmCell *cell = chip->cells->data[0]; - // readouts->n > 1 is not well-defined - if (cell->readouts->n > 1) { - psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nReadouts > 1"); - return false; + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!input) { + psErrorStackPrint(stderr, "Can't find input data!\n"); + exit(EXIT_FAILURE); } - pmReadout *readout = cell->readouts->data[0]; - psMetadata *header = pmReadoutGetHeader (readout); - // load the corresponding reference data (DVO command) - // XXX EAM : this needs to be improved, perhaps moved to a higher level - psArray *refstars = psastroLoadReference (header, config); + double plateScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLATE.SCALE"); + if (!status) plateScale = 1.0; - // we require PASTRO.OBJECTS (pmAstromObj) on the analysis for this readout - psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.OBJECTS"); + pmFPAview *view = pmFPAviewAlloc (0); + pmFPA *fpa = input->fpa; - // project the rawstars to the current best guess astrometry - psastroProjectRawstars (rawstars, readout); + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace (__func__, 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace (__func__, 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } - // use the header & config info to project refstars onto the focal plane - psastroProjectRefstars (refstars, readout); + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } - // testWriteRaw ("ref.inp", refstars); - // testWriteRaw ("raw.inp", rawstars); + // select the raw objects for this readout + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS"); + if (rawstars == NULL) { continue; } - // fprintf (stderr, "rawstars:\n"); - // psastroDumpStars (rawstars); - // fprintf (stderr, "refstars:\n"); - // psastroDumpStars (refstars); + // the refstars is a subset within range of this chip + psArray *refstars = psArrayAlloc (100); - // find initial offset / rotation - stats = pmAstromGridMatch (rawstars, refstars, config); + // select the reference objects within range of this readout + // project the reference objects to this chip + for (int i = 0; i < refs->n; i++) { + pmAstromObj *ref = refs->data[i]; + + p_psProject (ref->TP, ref->sky, fpa->projection); + psPlaneDistortApply (ref->FP, fpa->fromTangentPlane, ref->TP, 0.0, 0.0); + psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); - // adjust the chip.toFPA terms only - pmAstromGridApply (chip->toFPA, stats); - chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); + // XXX what is the X,Y range of the selected chip? + // XXX for the moment, force these to be good for Megacam + if (ref->chip->x < -500) continue; + if (ref->chip->x > 2500) continue; + if (ref->chip->y < -500) continue; + if (ref->chip->y > 5000) continue; + + psArrayAdd (refstars, 100, ref); + } + psastroOneChip (fpa, chip, refstars, rawstars, recipe); + psFree (refstars); - // use fit result to re-project rawstars - psastroProjectRawstars (rawstars, readout); - psastroProjectRefstars (refstars, readout); - - // testWriteRaw ("ref.dat", refstars); - // testWriteRaw ("raw.dat", rawstars); - - // use small radius to match stars - match = pmAstromRadiusMatch (rawstars, refstars, config); - - // improved fit for astrometric terms - pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, config); - chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); - - // per-chip astrometry has been updated - pmAstromWriteWCS (chip->toFPA, fpa->toSky, header); - + // read WCS data from the corresponding header + pmHDU *hdu = pmFPAviewThisHDU (view, fpa); + pmAstromWriteWCS (chip->toFPA, fpa->toSky, hdu->header, plateScale); + } + } + } + psFree (view); return true; } + +/* coordinate frame hierachy + pixels (on a given readout) + cell + chip + FP (focal plane) + TP (tangent plane) + sky (ra, dec) +*/ Index: /trunk/psastro/src/psastroConvert.c =================================================================== --- /trunk/psastro/src/psastroConvert.c (revision 7013) +++ /trunk/psastro/src/psastroConvert.c (revision 7014) @@ -1,4 +1,21 @@ # include "psastro.h" +// XXX leak free 2006.04.27 +// XXX pass/apply the WCS information? +bool psastroConvert (pmReadout *readout, psMetadata *recipe) { + + psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); + if (sources == NULL) + return false; + + psArray *objects = pmSourceToAstromObj (sources); + + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.OBJECTS", PS_DATA_ARRAY, "astrometry objects", objects); + psFree (objects); + + return true; +} + +// XXX select a magnitude range? psArray *pmSourceToAstromObj (psArray *sources) { @@ -6,8 +23,8 @@ for (int i = 0; i < sources->n; i++) { - psSource *source = sources->data[i]; + pmSource *source = sources->data[i]; // only accept the PSF sources? - if (source->type != PM_SOURCE_STAR) continue; + if (source->type != PM_SOURCE_TYPE_STAR) continue; pmModel *model = source->modelPSF; @@ -25,9 +42,10 @@ // XXX do we have the information giving the readout and cell offset? // for the moment, assume chip == cell == readout - obj->cell = obj->pix; - obj->chip = obj->cell; - - // can we make an initial guess of the astrometry to set FP, TP, and sky coords? - - + *obj->cell = *obj->pix; + *obj->chip = *obj->cell; + + psArrayAdd (objects, 100, obj); + psFree (obj); + } + return objects; } Index: /trunk/psastro/src/psastroDataLoad.c =================================================================== --- /trunk/psastro/src/psastroDataLoad.c (revision 7014) +++ /trunk/psastro/src/psastroDataLoad.c (revision 7014) @@ -0,0 +1,65 @@ +# include "psastro.h" +// XXX leak free 2006.04.27 + +// this loop loads the data from the input files and selects the +// brighter stars for astrometry +// at the end of this function, the complete stellar data is loaded +// into the correct fpa structure locations (readout.analysis:PSPHOT.SOURCES) +// +// all of the different astrometry analysis modes use the same data load loop +bool psastroDataLoad (pmConfig *config) { + + pmChip *chip; + pmCell *cell; + pmReadout *readout; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO"); + if (!recipe) { + psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n"); + exit(EXIT_FAILURE); + } + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!input) { + psErrorStackPrint(stderr, "Can't find input data!\n"); + exit(EXIT_FAILURE); + } + // de-activate all files except PSASTRO.INPUT + pmFPAfileActivate (config->files, false, NULL); + pmFPAfileActivate (config->files, true, "PSASTRO.INPUT"); + + pmFPAview *view = pmFPAviewAlloc (0); + + // files associated with the science image + pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); + + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psTrace (__func__, 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); + + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psTrace (__func__, 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); + if (! readout->data_exists) { continue; } + + psastroConvert (readout, recipe); + + pmFPAfileIOChecks (config->files, view, PM_FPA_AFTER); + } + pmFPAfileIOChecks (config->files, view, PM_FPA_AFTER); + } + pmFPAfileIOChecks (config->files, view, PM_FPA_AFTER); + } + pmFPAfileIOChecks (config->files, view, PM_FPA_AFTER); + psFree (view); + return true; +} + Index: unk/psastro/src/psastroDataLoop.c =================================================================== --- /trunk/psastro/src/psastroDataLoop.c (revision 7013) +++ (revision ) @@ -1,49 +1,0 @@ -# include "psastro.h" - -// this loop loads the data from the input files and selects the -// brighter stars for astrometry -// at the end of this function, the complete stellar data is loaded -// into the correct fpa structure locations (readout.analysis:PSPHOT.SOURCES) -// -// all of the different astrometry analysis modes use the same data load loop -bool psastroDataLoop (pmConfig *config) { - - pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSASTRO.INPUT"); - if (!status) { - psErrorStackPrint(stderr, "Can't find input data!\n"); - exit(EXIT_FAILURE); - } - // de-activate all files except PSASTRO.INPUT - psFPAfileActivate (config->files, false, NULL); - psFPAfileActivate (config->files, true, "PSASTRO.INPUT"); - - pmFPAview *view = pmFPAviewAlloc (0); - - // files associated with the science image - pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); - - while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); - if (!chip->process || !chip->file_exists) { continue; } - pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); - - while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { - psLogMsg ("psphot", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); - if (!cell->process || !cell->file_exists) { continue; } - pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); - - // process each of the readouts - while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { - pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); - if (! readout->data_exists) { continue; } - - pmFPAfileIOChecks (config->files, view, PM_FPA_AFTER); - } - pmFPAfileIOChecks (config->files, view, PM_FPA_AFTER); - } - pmFPAfileIOChecks (config->files, view, PM_FPA_AFTER); - } - pmFPAfileIOChecks (config->files, view, PM_FPA_AFTER); - return true; -} - Index: /trunk/psastro/src/psastroDataSave.c =================================================================== --- /trunk/psastro/src/psastroDataSave.c (revision 7013) +++ /trunk/psastro/src/psastroDataSave.c (revision 7014) @@ -1,40 +1,47 @@ # include "psastro.h" +// XXX leak free 2006.04.27 -// load the data from the files in this loop. -// we write out the result in a second loop -// at the end of this function, the complete stellar data is loaded -// into the correct fpa structure locations (readout.analysis:PSPHOT.SOURCES) +// this loop saves the photometry/astrometry data files bool psastroDataSave (pmConfig *config) { - // define the output files: - pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSASTRO.INPUT"); - if (!status) { - psErrorStackPrint(stderr, "Can't find input data!\n"); + pmChip *chip; + pmCell *cell; + pmReadout *readout; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO"); + if (!recipe) { + psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n"); exit(EXIT_FAILURE); } - // this is the output target - pmFPAfileDefine (config->files, format, input, "PSASTRO.OUTPUT"); + // select the output data sources + pmFPAfile *output = psMetadataLookupPtr (NULL, config->files, "PSASTRO.OUTPUT"); + if (!output) { + psErrorStackPrint(stderr, "Can't find output data!\n"); + exit(EXIT_FAILURE); + } - char *output = psMetadataLookupPtr(&status, config->arguments, "OUTPUT"); - pmFPAfileAddOutputName (config->files, "OUTPUT", output, PM_FPA_MODE_WRITE); + // de-activate all files except PSASTRO.OUTPUT + pmFPAfileActivate (config->files, false, NULL); + pmFPAfileActivate (config->files, true, "PSASTRO.OUTPUT"); pmFPAview *view = pmFPAviewAlloc (0); - // files associated with the science image + // open/load files as needed pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); - while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + while ((chip = pmFPAviewNextChip (view, output->fpa, 1)) != NULL) { + psTrace (__func__, 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); if (!chip->process || !chip->file_exists) { continue; } pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); - while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { - psLogMsg ("psphot", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + while ((cell = pmFPAviewNextCell (view, output->fpa, 1)) != NULL) { + psTrace (__func__, 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); if (!cell->process || !cell->file_exists) { continue; } pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); // process each of the readouts - while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + while ((readout = pmFPAviewNextReadout (view, output->fpa, 1)) != NULL) { pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE); if (! readout->data_exists) { continue; } @@ -47,4 +54,6 @@ } pmFPAfileIOChecks (config->files, view, PM_FPA_AFTER); + + psFree (view); return true; } Index: unk/psastro/src/psastroIO.c =================================================================== --- /trunk/psastro/src/psastroIO.c (revision 7013) +++ (revision ) @@ -1,230 +1,0 @@ -# include "psastro.h" - -// crude function to load CMP data (header + ascii) -// XXX EAM : assumes fixed line, expects NSTARS entries -psArray *psastroReadCMP (psMetadata **header, char *filename) { - - bool status; - char line[80]; - double X, Y, Mag, dMag; - psMetadata *myHeader; - - psFits *fits = psFitsAlloc (filename); - myHeader = psFitsReadHeader (NULL, fits); - psFree (fits); - - // how many lines in the header? - // XXX EAM : is this calculation robust? - int nLines = myHeader->list->n; - int nBytes = nLines * 80; - if (nBytes % 2880) { - int nBlock = 1 + (int)(nBytes / 2880); - nBytes = nBlock * 2880; - } - - // re-open, seek to end of header - FILE *f = fopen (filename, "r"); - if (f == NULL) { - psLogMsg ("pmSourcesReadCMP", 3, "can't open output file for input %s\n", filename); - return NULL; - } - fseek (f, nBytes, SEEK_SET); - - // prepare array to store data - int nStars = psMetadataLookupS32 (&status, myHeader, "NSTARS"); - psArray *stars = psArrayAlloc (nStars); - stars->n = 0; - - - for (int i = 0; i < nStars; i++) { - fread (line, 1, 66, f); - line[66] = 0; - sscanf (line, "%lf %lf %lf %lf", &X, &Y, &Mag, &dMag); - // fprintf (stderr, "%f %f %f %f\n", X, Y, Mag, dMag); - - pmAstromObj *obj = pmAstromObjAlloc (); - obj->pix.x = X; - obj->pix.y = Y; - obj->Mag = Mag; - obj->dMag = dMag; - psArrayAdd (stars, 100, obj); - } - fclose (f); - - *header = myHeader; - return (stars); -} - -// elixir-style pseudo FITS table (header + ascii list) -bool psastroWriteCMP (pmFPA *fpa, char *filename) { - - bool status; - - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - for (int j = 0; j < chip->cells->n; j++) { - pmCell *cell = chip->cells->data[j]; - - // create file, write-out header - unlink (filename); - psFits *fits = psFitsAlloc (filename); - psFitsWriteHeader (cell->header, fits); - psFree (fits); - - // re-open, add data to end of file - FILE *f = fopen (filename, "a+"); - if (f == NULL) { - psLogMsg ("psastroWriteCMP", 3, "can't open output file for output %s\n", filename); - return false; - } - fseek (f, 0, SEEK_END); - - for (int k = 0; k < cell->readouts->n; k++) { - pmReadout *readout = cell->readouts->data[k]; - psArray *stars = psMetadataLookupPtr (&status, readout->analysis, "STARS.FULLSET"); - for (int m = 0; m < stars->n; m++) { - pmAstromObj *star = stars->data[m]; - fprintf (f, "%6.1f %6.1f\n", star->pix.x, star->pix.y); - } - } - fclose (f); - } - } - return true; -} - -// fake this by loading from a text file ?? -psArray *psastroLoadReference (psMetadata *header, psMetadata *config) { - - // which reference source? - // ra & dec, area from header - - double RA, DEC, Mag, dMag; - - FILE *f = fopen ("test.cat", "r"); - if (f == NULL) { - psLogMsg ("psastroLoadReference", 3, "can't open dummy reference file test.ref\n"); - return NULL; - } - - psArray *stars = psArrayAlloc (100); - stars->n = 0; - - // we have fixed bytes / line : use that info - while (fscanf (f, "%lf %lf %lf %lf", &RA, &DEC, &Mag, &dMag) != EOF) { - pmAstromObj *obj = pmAstromObjAlloc (); - obj->sky.r = RA; - obj->sky.d = DEC; - obj->Mag = Mag; - obj->dMag = dMag; - psArrayAdd (stars, 100, obj); - } - fclose (f); - return (stars); -} - -// elixir-style pseudo FITS table (header + ascii list) -bool testWriteCMP (psMetadata *header, char *filename, psArray *sources) { - - int i; - psMetadataItem *mdi; - - psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "foo", sources->n); - - // create file, write-out header - unlink (filename); - psFits *fits = psFitsAlloc (filename); - - // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0 - mdi = psMetadataLookup (header, "NAXIS"); - mdi->data.S32 = 0; - mdi->type = PS_DATA_S32; - - psFitsWriteHeader (header, fits); - psFree (fits); - - // re-open, add data to end of file - FILE *f = fopen (filename, "a+"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); - return false; - } - fseek (f, 0, SEEK_END); - - for (i = 0; i < sources->n; i++) { - - pmAstromObj *star = sources->data[i]; - - fprintf (f, "%6.1f %6.1f %6.3f %03d %2d %3.1f %6.3f %6.3f %6.2f %6.2f %5.1f\n", - star->pix.x, star->pix.y, star->Mag, (int)(star->dMag*1000), 1, 2.5, 11.0, 12.0, 2.0, 0.5, 15.0); - } - fclose (f); - return true; -} - -// write out raw objects -bool testWriteRaw (char *filename, psArray *sources) { - - int i; - - // re-open, add data to end of file - FILE *f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); - return false; - } - - for (i = 0; i < sources->n; i++) { - - pmAstromObj *star = sources->data[i]; - - fprintf (f, "%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %10.6f %10.6f %8.2f %8.2f\n", - star->pix.x, star->pix.y, - star->cell.x, star->cell.y, - star->chip.x, star->chip.y, - star->FP.x, star->FP.y, - star->TP.x, star->TP.y, - star->sky.r*DEG_RAD, star->sky.d*DEG_RAD, - star->Mag, star->dMag); - } - fclose (f); - return true; -} - -// elixir-style pseudo FITS table (header + ascii list) -bool testWriteRef (char *filename, psArray *sources) { - - int i; - - // re-open, add data to end of file - FILE *f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); - return false; - } - - for (i = 0; i < sources->n; i++) { - - pmAstromObj *star = sources->data[i]; - - fprintf (f, "%10.7f %10.7f %6.3f %6.3f\n", star->sky.r, star->sky.d+10*RAD_DEG/3600, 10.0, 0.05); - } - fclose (f); - return true; -} - -void psastroDumpStars (psArray *sources) { - - for (int i = 0; i < sources->n; i++) { - pmAstromObj *star = sources->data[i]; - - fprintf (stderr, "%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %10.6f %10.6f %8.2f %8.2f\n", - star->pix.x, star->pix.y, - star->cell.x, star->cell.y, - star->chip.x, star->chip.y, - star->FP.x, star->FP.y, - star->TP.x, star->TP.y, - star->sky.r*DEG_RAD, star->sky.d*DEG_RAD, - star->Mag, star->dMag); - } -} Index: unk/psastro/src/psastroLoadData.c =================================================================== --- /trunk/psastro/src/psastroLoadData.c (revision 7013) +++ (revision ) @@ -1,38 +1,0 @@ -# include "psphot.h" - -bool psastroLoadData (ppFile *input, psDB *db, int chipNum, int cellNum) -{ - // an input chip is valid for processing if: - // (((chipNum == i) || (chipNum == -1)) && chip.process) - - // If we have not opened the file, skip it - if (input->fits == NULL) { - return false; - } - - psTrace(__func__, 1, "Loading %d,%d for %s\n", chipNum, cellNum, input->filename); - - // set input:valid flags according to process and chipNum/cellNum - for (int i = 0; i < input->fpa->chips->n; i++) { - pmChip *chip = input->fpa->chips->data[i]; // Chip in input image - chip->process = (! chip->exists && ((chipNum == i) || (chipNum == -1))); - - for (int j = 0; j < chip->cells->n; j++) { - pmCell *cell = chip->cells->data[j]; // Cell in input image - cell->process &= chip->process; - cell->process = (! cell->exists && ((cellNum == i) || (cellNum == -1))); - } - } - - // Read in the input pixels - if (! pmFPARead(input->fpa, input->fits, input->phu, db)) { - psErrorStackPrint(stderr, "Unable to populate camera from input FITS file\n"); - exit(EXIT_FAILURE); - } - - return true; -} - -// XXX this is not very efficient with fseeks : each pmFPARead is randomly accessing the file -// XXX does this handle multi-file data? -// XXX this does NOT preserve the state of the input valid flags Index: /trunk/psastro/src/psastroLoadReferences.c =================================================================== --- /trunk/psastro/src/psastroLoadReferences.c (revision 7014) +++ /trunk/psastro/src/psastroLoadReferences.c (revision 7014) @@ -0,0 +1,92 @@ +# include "psastro.h" +int mkstemp(char *template); + +# define ELIXIR_MODE 1 + +psArray *psastroLoadReferences (pmConfig *config) { + + int fd; + bool status; + char tempFile[64], tempLine[256]; + + // select the DVO database? + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO"); + + // DVO APIs expect decimal degrees + float RAmin = DEG_RAD*psMetadataLookupF32(NULL, recipe, "RA_MIN"); + float RAmax = DEG_RAD*psMetadataLookupF32(NULL, recipe, "RA_MAX"); + float DECmin = DEG_RAD*psMetadataLookupF32(NULL, recipe, "DEC_MIN"); + float DECmax = DEG_RAD*psMetadataLookupF32(NULL, recipe, "DEC_MAX"); + + // XXX need to add a padding area + // issue the following command: + // getstar -region RAmin RAmax DECmin DECmax + + sprintf (tempFile, "/tmp/psastro.XXXXXX"); + if ((fd = mkstemp (tempFile)) == -1) { + fprintf (stderr, "error creating temp output file\n"); + return NULL; + } + close (fd); + + // use fork to add timeout capability + // use cfitsio |filename format to avoid tempFile + # if ELIXIR_MODE + sprintf (tempLine, "getstar -D CATFORMAT elixir -D CATMODE mef -region %f %f %f %f > %s", RAmin, DECmin, RAmax, DECmax, tempFile); + # else + sprintf (tempLine, "getstar -D CATFORMAT panstarrs -D CATMODE mef -region %f %f %f %f > %s", RAmin, DECmin, RAmax, DECmax, tempFile); + # endif + + psTrace (__func__, 3, "%s\n", tempLine); + status = system (tempLine); + if (status) { + fprintf (stderr, "error loading reference data\n"); + return NULL; + } + + // the output from getstar is a file with the Average table + psFits *fits = psFitsOpen (tempFile, "r"); + + # if ELIXIR_MODE + psFitsMoveExtName (fits, "DVO_AVERAGE_ELIXIR"); + # else + psFitsMoveExtName (fits, "DVO_AVERAGE_PANSTARRS"); + # endif + + psMetadata *header = psFitsReadHeader (NULL, fits); + psArray *table = psFitsReadTable (fits); + psFitsClose (fits); + unlink (tempFile); + + // convert the Average table to the pmAstromObj entries + psArray *objects = psArrayAlloc (table->n); + for (int i = 0; i < table->n; i++) { + pmAstromObj *object = pmAstromObjAlloc (); + + psMetadata *row = table->data[i]; + + // DVO tables are stored in degrees + # if ELIXIR_MODE + object->sky->r = RAD_DEG*psMetadataLookupF32 (&status, row, "RA"); + object->sky->d = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC"); + object->Mag = 0.001*psMetadataLookupS32 (&status, row, "MAG"); // XXX ELIXIR uses millimags, PANSTARRS uses mags + # else + object->sky->r = RAD_DEG*psMetadataLookupF64 (&status, row, "RA"); + object->sky->d = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC"); + object->Mag = psMetadataLookupF32 (&status, row, "MAG"); + # endif + + psArrayAdd (objects, 100, object); + psFree (object); + } + psFree (header); + psFree (table); + + psTrace (__func__, 3, "loaded %d reference stars from (%10.6f,%10.6f) - (%10.6f,%10.6f)\n", + objects->n, RAmin, DECmin, RAmax, DECmax); + return objects; +} + +// XXX this function needs clarifiction of the getstar output format Index: /trunk/psastro/src/psastroOneChip.c =================================================================== --- /trunk/psastro/src/psastroOneChip.c (revision 7014) +++ /trunk/psastro/src/psastroOneChip.c (revision 7014) @@ -0,0 +1,39 @@ +# include "psastro.h" + +bool psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refstars, psArray *rawstars, psMetadata *recipe) { + + // bool status; + psArray *match; + + psastroWriteStars ("raw.0.dat", rawstars); + psastroWriteStars ("ref.0.dat", refstars); + psastroWriteTransform (chip->toFPA); + + // find initial offset / rotation + pmAstromStats stats = pmAstromGridMatch (rawstars, refstars, recipe); + stats = pmAstromGridTweak (rawstars, refstars, recipe, stats); + psLogMsg (__func__, 3, "grid search result: %f,%f @ %f deg\n", stats.offset.x, stats.offset.y, DEG_RAD*stats.angle); + + // adjust the chip.toFPA terms only + pmAstromGridApply (chip->toFPA, stats); + psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + + psastroWriteStars ("raw.1.dat", rawstars); + psastroWriteStars ("ref.1.dat", refstars); + psastroWriteTransform (chip->toFPA); + + // use small radius to match stars + match = pmAstromRadiusMatch (rawstars, refstars, recipe); + + // improved fit for astrometric terms + pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, recipe); + psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + + psastroWriteStars ("raw.2.dat", rawstars); + psastroWriteStars ("ref.2.dat", refstars); + psastroWriteTransform (chip->toFPA); + + psFree (match); + return true; +} + Index: /trunk/psastro/src/psastroParseCamera.c =================================================================== --- /trunk/psastro/src/psastroParseCamera.c (revision 7013) +++ /trunk/psastro/src/psastroParseCamera.c (revision 7014) @@ -1,10 +1,11 @@ -# include "psphot.h" +# include "psastro.h" +// XXX leak free 2006.04.27 -bool *psastroParseCamera (pmConfig *config) { +bool psastroParseCamera (pmConfig *config) { - ppFile *input = ppFileAlloc (); + bool status = false; // the input image(s) are required arguments; they define the camera - input = pmFPAfileAddSource (config, "INPUT", "PSASTRO.INPUT", true); + pmFPAfile *input = pmFPAfileFromArgs (&status, config, "PSASTRO.INPUT", "INPUT"); if (!status) { psAbort (__func__, "missing INPUT entry"); } @@ -15,6 +16,6 @@ // set default recipe values here - psMetadataAddStr (recipe, PS_LIST_TAIL, "PHOTCODE", PS_META_NO_REPLACE, "default fitting mode", "NONE"); - psMetadataAddStr (recipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_NO_REPLACE, "default fitting mode", "NONE"); + psMetadataAddStr (recipe, PS_LIST_TAIL, "PHOTCODE", PS_META_NO_REPLACE, "", "NONE"); + psMetadataAddStr (recipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_NO_REPLACE, "", "NONE"); // these calls bind the I/O handle to the specified fpa @@ -38,6 +39,22 @@ } } + psFree (chips); psTrace(__func__, 1, "Done with psastroParseCamera...\n"); return true; } + + +// useful for debugging +# if 0 +for (int i = 0; i < input->fpa->chips->n; i++) { + pmChip *chip = input->fpa->chips->data[i]; + fprintf (stderr, "chip %2d: %x %x\n", i, chip->file_exists, chip->process); + + for (int j = 0; j < chip->cells->n; j++) { + pmCell *cell = chip->cells->data[j]; + fprintf (stderr, "cell %2d: %x %x\n", j, cell->file_exists, cell->process); + + } +} +# endif Index: /trunk/psastro/src/psastroTestFuncs.c =================================================================== --- /trunk/psastro/src/psastroTestFuncs.c (revision 7014) +++ /trunk/psastro/src/psastroTestFuncs.c (revision 7014) @@ -0,0 +1,45 @@ +# include "psastro.h" + +// write out objects +bool psastroWriteStars (char *filename, psArray *sources) { + + // re-open, add data to end of file + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("psastroWriteStars", 3, "can't open output file for output %s\n", filename); + return false; + } + + for (int i = 0; i < sources->n; i++) { + + pmAstromObj *star = sources->data[i]; + + fprintf (f, "%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %10.6f %10.6f %8.2f %8.2f\n", + star->chip->x, star->chip->y, + star->FP->x, star->FP->y, + star->TP->x, star->TP->y, + star->sky->r*DEG_RAD, star->sky->d*DEG_RAD, + star->Mag, star->dMag); + } + fclose (f); + return true; +} + +bool psastroWriteTransform (psPlaneTransform *map) { + + // dump initial values: + for (int i = 0; i < map->x->nX + 1; i++) { + for (int j = 0; j < map->x->nY + 1; j++) { + if (map->x->mask[i][j]) continue; + psLogMsg (__func__, 4, "x term %d,%d: %f +/- %f\n", i, j, map->x->coeff[i][j], map->x->coeffErr[i][j]); + } + } + + for (int i = 0; i < map->y->nX + 1; i++) { + for (int j = 0; j < map->y->nY + 1; j++) { + if (map->y->mask[i][j]) continue; + psLogMsg (__func__, 4, "y term %d,%d: %f +/- %f\n", i, j, map->y->coeff[i][j], map->y->coeffErr[i][j]); + } + } + return true; +} Index: /trunk/psastro/src/psastroUtils.c =================================================================== --- /trunk/psastro/src/psastroUtils.c (revision 7013) +++ /trunk/psastro/src/psastroUtils.c (revision 7014) @@ -1,4 +1,29 @@ # include "psastro.h" # define RENORM 0 + +bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars) { + + // XXX this region needs to be defined more intelligently + psRegion region = {0, 0, 2000, 4500}; + psFree (chip->fromFPA); + chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, region, 20); + + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *raw = rawstars->data[i]; + + psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); + psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0); + p_psDeproject (raw->sky, raw->TP, fpa->projection); + } + + for (int i = 0; i < refstars->n; i++) { + pmAstromObj *ref = refstars->data[i]; + psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); + } + + return true; +} + +# if 0 bool psastroSelectBrightStars (pmFPA *fpa, psMetadata *config) { @@ -143,2 +168,146 @@ return true; } + +psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input) { + + psPolynomial2D *output = psPolynomial2DAlloc (input->nX, input->nY, input->type); + + for (int i = 0; i < input->nX; i++) { + for (int j = 0; j < input->nY; j++) { + output->mask[i][j] = input->mask[i][j]; + output->coeff[i][j] = input->coeff[i][j]; + output->coeffErr[i][j] = input->coeffErr[i][j]; + } + } + return (output); +} + +psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input) { + + psPolynomial4D *output = psPolynomial4DAlloc (input->nX, input->nY, input->nZ, input->nT, input->type); + + for (int i = 0; i < input->nX; i++) { + for (int j = 0; j < input->nY; j++) { + for (int k = 0; k < input->nZ; k++) { + for (int m = 0; m < input->nT; m++) { + output->mask[i][j][k][m] = input->mask[i][j][k][m]; + output->coeff[i][j][k][m] = input->coeff[i][j][k][m]; + output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m]; + } + } + } + } + return (output); +} + +psPlaneDistort *psPlaneDistortCopy (psPlaneDistort *input) { + + psPlaneDistort *output = psPlaneDistortAlloc (input->x->nX, input->x->nY, input->x->nZ, input->x->nT); + + for (int i = 0; i < input->x->nX; i++) { + for (int j = 0; j < input->x->nY; j++) { + for (int k = 0; k < input->x->nZ; k++) { + for (int m = 0; m < input->x->nT; m++) { + // x-terms + output->x->mask[i][j][k][m] = input->x->mask[i][j][k][m]; + output->x->coeff[i][j][k][m] = input->x->coeff[i][j][k][m]; + output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m]; + // y-terms + output->y->mask[i][j][k][m] = input->y->mask[i][j][k][m]; + output->y->coeff[i][j][k][m] = input->y->coeff[i][j][k][m]; + output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m]; + } + } + } + } + return (output); +} + +psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input) { + + psPlaneTransform *output = psPlaneTransformAlloc (input->x->nX, input->x->nY); + + for (int i = 0; i < input->x->nX; i++) { + for (int j = 0; j < input->x->nY; j++) { + // x-terms + output->x->mask[i][j] = input->x->mask[i][j]; + output->x->coeff[i][j] = input->x->coeff[i][j]; + output->x->coeffErr[i][j] = input->x->coeffErr[i][j]; + // y-terms + output->y->mask[i][j] = input->y->mask[i][j]; + output->y->coeff[i][j] = input->y->coeff[i][j]; + output->y->coeffErr[i][j] = input->y->coeffErr[i][j]; + } + } + return (output); +} + +psProjection *psProjectionCopy (psProjection *input) { + + psProjection *output = psProjectionAlloc (input->R, input->D, input->Xs, input->Ys, input->type); + return (output); +} + +// very crude distortion inversion: assumes 0 order in z and t, linear in x and y: +psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort) { + PS_ASSERT_PTR_NON_NULL(distort, 0); + PS_ASSERT_PTR_NON_NULL(distort->x, 0); + PS_ASSERT_PTR_NON_NULL(distort->y, 0); + + psPlaneDistort *out = psPlaneDistortAlloc(1, 1, 0, 0); + + /* simple matrix inversion code */ + + psF64 r11 = distort->x->coeff[1][0][0][0]; + psF64 r12 = distort->x->coeff[0][1][0][0]; + psF64 r21 = distort->y->coeff[1][0][0][0]; + psF64 r22 = distort->y->coeff[0][1][0][0]; + psF64 xo = distort->x->coeff[0][0][0][0]; + psF64 yo = distort->y->coeff[0][0][0][0]; + + psF64 invDet = 1.0 / (r11 * r22 - r12 * r21); // Inverse of the determinant + + out->x->coeff[1][0][0][0] = +invDet * r22; + out->x->coeff[0][1][0][0] = -invDet * r12; + out->y->coeff[1][0][0][0] = -invDet * r21; + out->y->coeff[0][1][0][0] = +invDet * r11; + + out->x->coeff[0][0][0][0] = - invDet * (r22 * xo - r12 * yo); + out->y->coeff[0][0][0][0] = - invDet * (r11 * yo - r21 * xo); + + return(out); +} + +// returns the rotation term, forcing positive parity +double psPlaneTransformGetRotation (psPlaneTransform *map) { + + if (map->x->nX < 1) return 0; + if (map->x->nY < 1) return 0; + + if (map->y->nX < 1) return 0; + if (map->y->nY < 1) return 0; + + double pc1_1 = map->x->coeff[1][0]; + double pc1_2 = map->x->coeff[0][1]; + double pc2_1 = map->y->coeff[1][0]; + double pc2_2 = map->y->coeff[0][1]; + + double px = SIGN (pc1_1); + double py = SIGN (pc2_2); + + // both x and y terms imply an angle. take the average + double t1 = -atan2 (px*pc1_2, px*pc1_1); + double t2 = +atan2 (py*pc2_1, py*pc2_2); + + // careful near -pi,+pi boundary... + if (t1 - t2 > M_PI/2) t2 += 2*M_PI; + if (t2 - t1 > M_PI/2) t1 += 2*M_PI; + + double theta = 0.5*(t1 + t2); + while (theta < M_PI) theta += 2*M_PI; + while (theta > M_PI) theta -= 2*M_PI; + + return (theta); +} + +# endif Index: /trunk/psastro/src/psastroWCS.c =================================================================== --- /trunk/psastro/src/psastroWCS.c (revision 7013) +++ /trunk/psastro/src/psastroWCS.c (revision 7014) @@ -2,7 +2,9 @@ // interpret header WCS (only handles traditional WCS for the moment) -bool pmAstromReadWCS (psPlaneTransform **toFPAOut, psProjection **toSkyOut, psMetadata *header) { +// plateScale is nominal physical scale on tangent plane (microns / arcsecond) +bool pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale) { psPlaneTransform *toFPA; + psPlaneDistort *toTPA; psProjection *toSky; psProjectionType type; @@ -19,4 +21,5 @@ // determine projection type + // XXX add the Elixir DIS / WRP two-layer projection here type = PS_PROJ_NTYPE; if (!strcmp (&ctype[4], "-SIN")) type = PS_PROJ_SIN; @@ -72,5 +75,5 @@ // renormalize to cdelt1, cdelt2, etc - float scale = hypot (pc1_1, pc1_2); + double scale = hypot (pc1_1, pc1_2); cdelt1 = cdelt2 = scale; pc1_1 /= scale; @@ -89,25 +92,44 @@ // scale = scale(i)/scale(0) (i == chip #) // project crval1(0),crval2(0 using projection - - // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581) - // XXX EAM : I've fixed this in pslib eam_rel8_b2 toFPA = psPlaneTransformAlloc (1, 1); - toFPA->x->coeff[0][0] = crpix1; - toFPA->x->coeff[1][0] = pc1_1; - toFPA->x->coeff[0][1] = pc1_2; + double cdelt = hypot (cdelt1, cdelt2) / plateScale; // degrees / micron (eg, in fact, whatever unit user chooses for focal plane) + cdelt1 /= cdelt; + cdelt2 /= cdelt; + + toFPA->x->coeff[0][0] = -(pc1_1*cdelt1*crpix1 + pc1_2*cdelt2*crpix2); + toFPA->x->coeff[1][0] = +(pc1_1*cdelt1); + toFPA->x->coeff[0][1] = +(pc1_2*cdelt2); toFPA->x->mask[1][1] = 1; - toFPA->y->coeff[0][0] = crpix2; - toFPA->y->coeff[1][0] = pc2_1; - toFPA->y->coeff[0][1] = pc2_2; + toFPA->y->coeff[0][0] = -(pc2_1*cdelt1*crpix1 + pc2_2*cdelt2*crpix2); + toFPA->y->coeff[1][0] = +(pc2_1*cdelt1); + toFPA->y->coeff[0][1] = +(pc2_2*cdelt2); toFPA->y->mask[1][1] = 1; - - // center of projection is (0,0) coordinate of TPA - toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type); - - *toFPAOut = toFPA; - *toSkyOut = toSky; - + *toFPAout = toFPA; + + if (*toSkyOut != NULL) { + if (*toTPAout == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined"); + + psSphere sky; + psPlane tpa; + + sky.r = crval1*RAD_DEG; + sky.d = crval2*RAD_DEG; + p_psProject (&tpa, &sky, *toSkyOut); + + // XXX for the moment, assume toTPA is the identity transformation + toFPA->x->coeff[0][0] = tpa.x; + toFPA->y->coeff[0][0] = tpa.y; + } else { + // XXX for now, use the identity for TPA <--> FPA + toTPA = psPlaneDistortIdentity (); + + // center of projection is (0,0) coordinate of TPA + toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, RAD_DEG*cdelt, RAD_DEG*cdelt, type); + + *toTPAout = toTPA; + *toSkyOut = toSky; + } return true; } @@ -115,22 +137,23 @@ // convert toFPA / toSky components to traditional WCS -bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header) { +// plateScale is nominal physical scale on tangent plane (microns / arcsecond) +bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale) { switch (toSky->type) { case PS_PROJ_SIN: - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---SIN"); - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--SIN"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---SIN"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--SIN"); break; case PS_PROJ_TAN: - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---TAN"); - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--TAN"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---TAN"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--TAN"); break; case PS_PROJ_AIT: - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---AIT"); - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--AIT"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---AIT"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--AIT"); break; case PS_PROJ_PAR: - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---PAR"); - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--PAR"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---PAR"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--PAR"); break; default: @@ -139,15 +162,18 @@ } - psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->R*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->D*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][0]); - psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][0]); - psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Xs*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Ys*DEG_RAD); - - psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[1][0]); - psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][1]); - psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[1][0]); - psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][1]); + // XXX not really right: needs to deal with non-identity toTP coeffs + // XXX actually, totally wrong. fix the conversions + // XXX need to handle the plateScale + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", toFPA->x->coeff[0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", toFPA->y->coeff[0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", toSky->Xs*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*DEG_RAD); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]); // alternative representations use @@ -157,2 +183,20 @@ return true; } + +psPlaneDistort *psPlaneDistortIdentity () { + + psPlaneDistort *distort; + + distort = psPlaneDistortAlloc (1, 1, 0, 0); + distort->x->coeff[0][0][0][0] = 0; + distort->x->coeff[1][0][0][0] = 1; + distort->x->coeff[0][1][0][0] = 0; + distort->x->mask [1][1][0][0] = 1; + + distort->y->coeff[0][0][0][0] = 0; + distort->y->coeff[1][0][0][0] = 0; + distort->y->coeff[0][1][0][0] = 1; + distort->y->mask [1][1][0][0] = 1; + + return distort; +} Index: unk/psastro/src/testPSastroLoop.c =================================================================== --- /trunk/psastro/src/testPSastroLoop.c (revision 7013) +++ (revision ) @@ -1,36 +1,0 @@ - -bool ImageLoop (psphotData *data, ppConfig *config) { - - psRegion region = {0,0,0,0}; - pmFPA *fpa = data->input->fpa; - pmFPAiterator *fpi = pmFPAiteratorAlloc (fpa, region); - - // the depths are probably interpreted from the camera config file - pmFileAlloc (fpi, config->camera, "PSASTRO.INPUT"); - pmFileAlloc (fpi, config->camera, "PSASTRO.OUTPUT"); - pmFileIOChecks (fpi); - - while ((chip = pmChipNext (fpi)) != NULL) { - - psLogMsg ("psphot", 4, "Chip %d: %x %x\n", i, chip->exists, chip->process); - if (! chip->process) { continue; } - pmFileIOChecks (fpi); - - while ((cell = pmCellNext (fpi)) != NULL) { - - psLogMsg ("psphot", 4, "Cell %d: %x %x\n", j, cell->exists, cell->process); - if (! cell->process) { continue; } - pmFileIOChecks (fpi); - - // process each of the readouts - while ((readout = pmReadoutNext (fpi) != NULL) { - pmFileIOChecks (fpi); - - // run the actual photometry analysis - psphotReadout (readout, fpi, config); - } - } - } - psphotOutputClose (data); - return true; -}