Index: /trunk/psastro/src/psModUtils.c =================================================================== --- /trunk/psastro/src/psModUtils.c (revision 5506) +++ /trunk/psastro/src/psModUtils.c (revision 5506) @@ -0,0 +1,115 @@ +# include "psastro.h" + +// these functions temporarily replace existing psModules code + +// template sample alloc/free pair: +static void psFooFree (psFoo *foo) { + + if (foo == NULL) return; + return; +} + +psFoo *psFooAlloc (int i1, int i2) { + + psFoo *foo = psAlloc (sizeof(psFoo)); + psMemSetDeallocator(foo, (psFreeFunc) psFooFree); + + return (foo); +} + +// pmFPA +static void pmFPAFree (pmFPA *fpa) { + + if (fpa == NULL) return; + + psFree (fpa->toSky); + psFree (fpa->toTPA); + psFree (fpa->fromTPA); + psFree (fpa->chips); + + return; +} + +pmFPA *pmFPAAlloc (int i1, int i2) { + + pmFPA *fpa = psAlloc (sizeof(pmFPA)); + psMemSetDeallocator(fpa, (psFreeFunc) pmFPAFree); + + fpa->toSky = NULL; + fpa->toTPA = NULL; + fpa->fromTPA = NULL; + fpa->chips = NULL; + + return (fpa); +} + +// pmChip +static void pmChipFree (pmChip *chip) { + + if (chip == NULL) return; + + psFree (chip->toFPA); + psFree (chip->fromFPA); + psFree (chip->cells); + + return; +} + +pmChip *pmChipAlloc (int i1, int i2) { + + pmChip *chip = psAlloc (sizeof(pmChip)); + psMemSetDeallocator(chip, (psFreeFunc) pmChipFree); + + chip->toFPA = NULL; + chip->fromFPA = NULL; + chip->cells = NULL; + + return (chip); +} + +// pmCell +static void pmCellFree (pmCell *cell) { + + if (cell == NULL) return; + + psFree (cell->toChip); + psFree (cell->readouts); + + return; +} + +pmCell *pmCellAlloc (int i1, int i2) { + + pmCell *cell = psAlloc (sizeof(pmCell)); + psMemSetDeallocator(cell, (psFreeFunc) pmCellFree); + + cell->toChip = NULL; + cell->readouts = NULL; + + return (cell); +} + +// pmReadout +static void pmReadoutFree (pmReadout *readout) { + + if (readout == NULL) return; + + psFree (readout->objects); + + return; +} + +pmReadout *pmReadoutAlloc (int i1, int i2) { + + pmReadout *readout = psAlloc (sizeof(pmReadout)); + psMemSetDeallocator(readout, (psFreeFunc) pmReadoutFree); + + readout->col0 = 0; + readout->row0 = 0; + readout->colBins = 1; + readout->rowBins = 1; + readout->objects = NULL; + + return (readout); +} + Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 5505) +++ /trunk/psastro/src/psastro.h (revision 5506) @@ -3,53 +3,47 @@ # include // for unlink # include -# include "psLibUtils.h" -# include "pmObjects_EAM.h" -# include "psModulesUtils.h" -# include "pmPSF.h" -# include "pmPSFtry.h" -# include "pmModelGroup.h" -typedef struct { - psImage *image; - psImage *mask; - psImage *weight; - psMetadata *header; -} eamReadout; +typedef struct +{ + psProjection *toSky; ///< Projection from tangent plane to sky + psPlaneDistort *toTPA; ///< Transformation from focal plane to tangent plane + psPlaneDistort *fromTPA; ///< Transformation from focal plane to tangent plane + psArray *chips; ///< The chips +} +pmFPA; -// top-level psphot functions -psMetadata *psphotArguments (int *argc, char **argv); -eamReadout *psphotSetup (psMetadata *config); -psStats *psphotImageStats (eamReadout *imdata, psMetadata *config); -psArray *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *allpeaks); -pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky); -bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats); -void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +typedef struct +{ + // Astrometric transformations + psPlaneTransform *toFPA; ///< Transformation from chip to FPA coordinates + psPlaneTransform *fromFPA; ///< Transformation from FPA to chip coordinates + psArray *cells; ///< The cells +} +pmChip; -bool psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE); -bool psphotSubtractPSF (pmSource *source); -int psphotSortBySN (const void **a, const void **b); -int psphotSaveImage (psMetadata *header, psImage *image, char *filename); -bool psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius); +typedef struct +{ + // Astrometric transformations + psPlaneTransform *toChip; ///< Transformations from cell to chip coordinates + psArray *readouts; ///< The readouts (referred to by number) + psMetadata *header; ///< header always corresponds to a cell +} +pmCell; -psArray *pmPeaksSigmaLimit (eamReadout *imdata, psMetadata *config, psStats *sky); -pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); -pmModel *pmSourceSelectModel (pmSource *source); +typedef struct +{ + // Position on the cell + int col0; ///< Offset from the left of chip. + int row0; ///< Offset from the bottom of chip. + int colBins; ///< Amount of binning in x-dimension + int rowBins; ///< Amount of binning in y-dimension + psArray *objects; ///< sources detected / measured on readout +} +pmReadout; -// eamReadout functions -eamReadout *eamReadoutAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header); - -// output functions -bool pmSourcesWriteText (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteOBJ (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteCMP (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteCMF (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); - -bool pmPeaksWriteText (psArray *sources, char *filename); -bool pmMomentsWriteText (psArray *sources, char *filename); -bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf); -bool pmModelWriteFLTs (psArray *sources, char *filename); -bool pmModelWriteNULLs (psArray *sources, char *filename); - -int pmSourcesDophotType (pmSource *source); +// dummy structure for template code (psFooAlloc, etc) +typedef struct +{ + int i; +} +psFoo; Index: /trunk/psastro/src/psastroBuildFPA.c =================================================================== --- /trunk/psastro/src/psastroBuildFPA.c (revision 5505) +++ /trunk/psastro/src/psastroBuildFPA.c (revision 5506) @@ -1,2 +1,175 @@ # include "psastro.h" +pmFPA *psastroBuildFPA (psMetadata *header, psMetadata *config) { + + // 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 (); + + fpa->chips = psArrayAlloc (Nchips); + for (int i = 0; i < Nchips; i++) { + + pmChip *chip = pmChipAlloc (); + chip->fpa = fpa; // assign parent fpa (view only; don't free) + + chip->cells = psArrayAlloc (Ncells); + for (int j = 0; j < Ncells; j++) { + + pmCell *cell = pmCellAlloc (); + cell->chip = chip; // assign parent chip (view only; don't free) + cell->header = header; + + pmCellInterpretWCS (cell, header); + + cell->readouts = psArrayAlloc (Nreadouts); + for (int k = 0; k < Nreadouts; k++) { + + pmReadout *readout = pmReadoutAlloc (); + cells->readouts->data[j] = readout; + } + chips->cells->data[j] = cell; + } + fpa->chips->data[i] = chip; + } + return (fpa); +} + +// interpret header WCS +bool pmCellInterpretWCS (pmCell *cell, psMetadata *header) { + + float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2; + float pc1_1, pc1_2, pc2_1, pc2_2; + + // *** interpret header data, convert to crval(i), etc + char *ctype = pmMetadataLookupPtr (&status, header, "CTYPE2"); + if (!status) { + psLogMsg ("psastro", 2, "warning: no WCS metadata in header\n"); + return false; + } + + // determine projection type + type = PS_PROJ_NTYPE; + if (!strcmp (&ctype[4], "-SIN")) type = PS_PROJ_SIN; + if (!strcmp (&ctype[4], "-TAN")) type = PS_PROJ_TAN; + if (!strcmp (&ctype[4], "-AIT")) type = PS_PROJ_AIT; + if (!strcmp (&ctype[4], "-PAR")) type = PS_PROJ_PAR; + if (type == PS_PROJ_NTYPE) { + psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype); + return false; + } + + crval1 = psMetadataLookupF32 (&status, header, "CRVAL1"); + crval2 = psMetadataLookupF32 (&status, header, "CRVAL2"); + crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1"); + crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2"); + + // test the CDELTi varient + cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1"); + if (status) { + cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2"); + + // test the CROTAi varient: + float rotate = psMetadataLookupF32 (&status, header, "CROTA2"); + if (status) { + Lambda = cdelt2 / cdelt1; + pc1_1 = cos(rotate*RAD_DEG); + pc1_2 = -sin(rotate*RAD_DEG) * Lambda; + pc2_1 = sin(rotate*RAD_DEG) / Lambda; + pc2_2 = cos(rotate*RAD_DEG); + goto got_matrix; + } + + // test the PC00i00j varient: + pc1_1 = psMetadataLookupF32 (&status, header, "PC001001"); + if (status) { + pc1_2 = psMetadataLookupF32 (&status, header, "PC001002"); + pc2_1 = psMetadataLookupF32 (&status, header, "PC002001"); + pc2_2 = psMetadataLookupF32 (&status, header, "PC002002"); + + // XXX EAM : add Elixir polynomial terms here eventually + goto got_matrix; + } + psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n"); + return false; + } + + // test the CDi_j varient + pc1_1 = psMetadataLookupF32 (&status, header, "CD1_1"); + if (status) { + pc1_2 = psMetadataLookupF32 (&status, header, "CD1_2"); + pc2_1 = psMetadataLookupF32 (&status, header, "CD2_1"); + pc2_2 = psMetadataLookupF32 (&status, header, "CD2_2"); + + // renormalize to cdelt1, cdelt2, etc + float scale = hypot (pc1_1, pc1_2); + cdelt1 = cdelt2 = scale; + pc1_1 /= scale; + pc1_2 /= scale; + pc2_1 /= scale; + pc2_2 /= scale; + goto got_matrix; + } + psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n"); + return false; + +got_matrix: + + // XXX EAM : these must already be set + chip = cell->chip; + fpa = chip->fpa; + + // set toChip to identity as default + // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581) + cell->toChip = psPlaneTransformAlloc (2, 2); + 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; + + // XXX EAM : if fpa->toSky and fpa->toTPA are already defined, then the + // toFPA must be modified to match the crval(i), scale(i) and crpix(i) + + // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581) + chip->toFPA = psPlaneTransformAlloc (2, 2); + + chip->toFPA->x->coeff[0][0] = crpix1; + chip->toFPA->x->coeff[1][0] = pc1_1; + chip->toFPA->x->coeff[0][1] = pc1_2; + chip->toFPA->x->mask[1][1] = 1; + + chip->toFPA->y->coeff[0][0] = crpix2; + chip->toFPA->y->coeff[1][0] = pc2_1; + chip->toFPA->y->coeff[0][1] = pc2_2; + chip->toFPA->y->mask[1][1] = 1; + + // set toTPA to identity as default + // XXX EAM : psPlaneDistortAlloc uses nTerm not nOrder (bug 581) + if (fpa->toTPA == NULL) { + fpa->toTPA = psPlaneDistortAlloc (2, 2, 1, 1); + 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; + } else { + psLogMsg ("psastro", 2, "warning: fpa distortion already defined\n"); + } + + // center of projection is (0,0) coordinate of TPA + fpa->toSky = psProjectionAlloc (crval1, crval2, cdelt1, cdelt2, type); + return true; +}