Changeset 5506

Timestamp:

Nov 11, 2005, 9:54:59 AM (20 years ago)

Author:

eugene

Message:

adding WCS interpretation code

Location:

trunk/psastro/src

Files:

• psModUtils.c (added)
• psastro.h (modified) (1 diff)
• psastroBuildFPA.c (modified) (1 diff)

trunk/psastro/src/psastro.h

@@ -3 +3 @@
 # include <unistd.h>   // for unlink
 # include <pslib.h>
-# 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;

trunk/psastro/src/psastroBuildFPA.c

@@ -1 +1 @@
 # 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;
+}