IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 7014


Ignore:
Timestamp:
Apr 30, 2006, 12:15:03 PM (20 years ago)
Author:
eugene
Message:

first complete working version (chip only)

Location:
trunk/psastro
Files:
5 added
12 deleted
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/Makefile

    r5620 r7014  
    2222PSASTRO = \
    2323$(SRC)/psastro.$(ARCH).o           \
    24 $(SRC)/psastroIO.$(ARCH).o         \
    25 $(SRC)/psastroBuildFPA.$(ARCH).o   \
    2624$(SRC)/psastroArguments.$(ARCH).o  \
    27 $(SRC)/psastroUtils.$(ARCH).o      \
    28 $(SRC)/psModUtils.$(ARCH).o        \
    29 $(SRC)/psLibUtils.$(ARCH).o        \
    30 $(SRC)/pmAstrom.$(ARCH).o          \
    31 $(SRC)/pmAstromGrid.$(ARCH).o
    32 
    33 PSASTRO-MKTEST = \
    34 $(SRC)/psastro-mktest.$(ARCH).o    \
    35 $(SRC)/psastroArguments.$(ARCH).o  \
    36 $(SRC)/psastroIO.$(ARCH).o         \
    37 $(SRC)/psModUtils.$(ARCH).o        \
    38 $(SRC)/psLibUtils.$(ARCH).o        \
    39 $(SRC)/pmAstrom.$(ARCH).o
     25$(SRC)/psastroParseCamera.$(ARCH).o\
     26$(SRC)/psastroDataLoad.$(ARCH).o   \
     27$(SRC)/psastroDataSave.$(ARCH).o   \
     28$(SRC)/psastroAstromGuess.$(ARCH).o\
     29$(SRC)/psastroLoadReferences.$(ARCH).o\
     30$(SRC)/psastroConvert.$(ARCH).o    \
     31$(SRC)/psastroChipAstrom.$(ARCH).o\
     32$(SRC)/psastroOneChip.$(ARCH).o\
     33$(SRC)/psastroUtils.$(ARCH).o\
     34$(SRC)/psastroTestFuncs.$(ARCH).o\
     35$(SRC)/psastroWCS.$(ARCH).o       
    4036
    4137psastro: $(BIN)/psastro.$(ARCH)
    4238$(BIN)/psastro.$(ARCH) : $(PSASTRO)
    43 $(PSASTRO): $(SRC)/psastro.h $(SRC)/pmAstrom.h
     39$(PSASTRO): $(SRC)/psastro.h
    4440
    45 psastro-mktest: $(BIN)/psastro-mktest.$(ARCH)
    46 $(BIN)/psastro-mktest.$(ARCH) : $(PSASTRO-MKTEST)
    47 $(PSASTRO-MKTEST): $(SRC)/psastro.h $(SRC)/pmAstrom.h
     41INSTALL = psastro
    4842
    49 INSTALL = psastro psastro-mktest
    50 
    51 all: psastro.install psastro-mktest.install
     43all: psastro.install
    5244
    5345# dependancy rules for binary code #########################

  • trunk/psastro/src/psastro.c

    r6792 r7014  
    55    psTimerStart ("complete");
    66
    7     psMetadata *header = NULL;
     7    // model inits are needed in pmSourceIO
     8    // XXX do we want to do the local model inits a la psphot?
     9    pmModelGroupInit ();
    810
    911    // load configuration information
    1012    pmConfig *config = psastroArguments (&argc, argv);
    1113
    12     // load input data (config and images (signal, noise, mask)
     14    // load identify the data sources
    1315    psastroParseCamera (config);
    1416
    15     // perform the astrometric solution
    16     psastroDataLoop (input, config);
     17    // load the raw pixel data (from PSPHOT.SOURCES)
     18    // select subset of stars for astrometry
     19    psastroDataLoad (config);
    1720
    18     // select astrometry stars from the psphot sources
    19     // limit fit to bright stars only
    20     psastroSelectAstromStars (config);
     21    // interpret the available initial astrometric information
     22    // apply the initial guess
     23    psastroAstromGuess (config);
    2124
    2225    // load the reference stars overlapping the data stars
    23     psastroLoadReference (config);
     26    psArray *refs = psastroLoadReferences (config);
    2427
    25     // fpa and subset point to the same astrometry terms
     28    psastroChipAstrom (config, refs);
     29
     30    # if 0
     31    // perform the desired astrometry analysis
    2632    switch (mode) {
    2733      case stack:
    28         psastroStackAstrom (config);
     34        psastroStackAstrom (config, refs);
     35        break
     36      case chip:
     37        psastroChipAstrom (config, refs);
    2938        break
    3039      case mosaic:
    31         psastroMosaicAstrom (config);
    32         break
    33       case mosaic:
    34         psastroMosaicAstrom (config);
     40        psastroMosaicAstrom (config, refs);
    3541        break;
    3642      default:
    3743        break;
    3844    }
     45    # endif
    3946
    40     // perform the astrometric solution
     47    // write out the results
    4148    psastroDataSave (config);
    4249
    43     psLogMsg ("psphot", 3, "complete psastro run: %f sec\n", psTimerMark ("complete"));
     50    psLogMsg ("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete"));
    4451
     52    psFree (refs);
    4553    psFree (config);
    46     psFree (fpa);
     54
     55    psTimerStop ();
     56    psMemCheckCorruption (true);
     57    pmModelGroupCleanup ();
     58    psTimeFinalize ();
     59    pmConceptsDone ();
     60    pmConfigDone ();
     61    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro");
     62    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro");
    4763
    4864    exit (0);

  • trunk/psastro/src/psastro.h

    r6791 r7014  
    33# include <unistd.h>   // for unlink
    44# include <pslib.h>
    5 # include <pmAstrometry.h>
    6 # include "psLibUtils.h"
    7 # include "pmAstrom.h"
     5# include <psmodules.h>
     6
     7# define psMemCopy(A)(psMemIncrRefCounter((A)))
     8# define DEG_RAD 57.295779513082322
     9# define RAD_DEG  0.017453292519943
     10# define SIGN(X)  (((X) == 0) ? 0 : ((fabs((double)(X))) / (X)))
    811
    912# define fromTPA fromTangentPlane
    … …  
    1114# define toSky projection
    1215
    13 // How much of the FPA to load at a time
    14 typedef enum {
    15     PP_LOAD_NONE,                       // Don't load anything
    16     PP_LOAD_FPA,                        // Load the entire FPA at once
    17     PP_LOAD_CHIP,                       // Load by chip
    18     PP_LOAD_CELL,                       // Load by cell
    19 } ppImageLoadDepth;
     16pmConfig         *psastroArguments (int *argc, char **argv);
     17bool              psastroParseCamera (pmConfig *config);
     18bool              psastroDataLoad (pmConfig *config);
     19bool              psastroDataSave (pmConfig *config);
    2020
    21 // Configuration data
    22 typedef struct {
    23     psMetadata *site;                   // The site configuration
    24     psMetadata *camera;                 // The camera configuration
    25     psMetadata *recipe;                 // The recipe (i.e., specific setups)
    26     psMetadata *arguments;              // Command-line arguments
    27     psMetadata *options;                // Command-line recipe options
    28     psDB       *database;               // Database handle
    29 } ppConfig;
     21bool              psastroConvert (pmReadout *readout, psMetadata *recipe);
     22psArray          *pmSourceToAstromObj (psArray *sources);
     23bool              psastroAstromGuess (pmConfig *config);
    3024
    31 // A file to process
    32 typedef struct {
    33     char *filename;                     // File name
    34     psFits *fits;                       // The FITS file handle
    35     psMetadata *phu;                    // The FITS header
    36     pmFPA *fpa;                         // The FPA, with pixels and extensions
    37 } ppFile;
     25bool              pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale);
     26bool              pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale);
     27psPlaneDistort   *psPlaneDistortIdentity ();
    3828
    39 bool              pmAstromReadWCS (psPlaneTransform **toFPA, psProjection **toSky, psMetadata *header);
    40 bool              pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header);
     29psArray          *psastroLoadReferences (pmConfig *config);
     30bool              psastroChipAstrom (pmConfig *config, psArray *refs);
     31bool              psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refset, psArray *rawset, psMetadata *recipe);
    4132
    42 pmFPA            *pmFPACopyAstrom (pmFPA *inFPA);
    43 psMetadata       *psastroArguments (int *argc, char **argv);
    44 pmFPA            *psastroBuildFPA (psMetadata *header, psArray *stars);
    45 bool              psastroChipAstrom (pmFPA *fpa, psMetadata *config);
    46 psArray          *psastroLoadReference (psMetadata *header, psMetadata *config);
    47 psArray          *psastroReadCMP (psMetadata **header, char *filename);
    48 bool              psastroSelectBrightStars (pmFPA *fpa, psMetadata *config);
    49 bool              psastroWriteCMP (pmFPA *fpa, char *filename);
     33// bool              psVectorSmooth (psVector *vector, double sigma, double Nsigma);
    5034
     35// utility functions:
     36bool              psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars);
     37bool              psastroWriteStars (char *filename, psArray *sources);
     38bool              psastroWriteTransform (psPlaneTransform *map);
    5139
    52 psPlane          *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout);
    53 psPlane          *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout);
    54 psPlane          *psCoordChipToCell_EAM(psPlane* cellCoord, const psPlane* chipCoord, const pmCell* cell);
    55 
    56 bool              testWriteCMP (psMetadata *header, char *filename, psArray *sources);
    57 bool              testWriteRef (char *filename, psArray *sources);
    58 bool              testWriteRaw (char *filename, psArray *sources);
    59 psMetadata       *testArguments (int *argc, char **argv);
    60 
    61 bool              psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky);
    62 bool              psastroProjectRawstars (psArray *stars, pmReadout *readout);
    63 bool              psastroProjectRefstars (psArray *stars, pmReadout *readout);
    64 
    65 psPlaneTransform *p_psPlaneTransformLinearInvert_EAM(psPlaneTransform *transform);
    66 pmFPA            *pmFPACopyAstrom (pmFPA *inFPA);

  • trunk/psastro/src/psastroArguments.c

    r6911 r7014  
    11# include "psastro.h"
    22# include <glob.h>
     3// XXX leak free 2006.04.27
    34
    45static void usage (void) {
    … …  
    910pmConfig *psastroArguments (int *argc, char **argv) {
    1011
     12    bool status;
    1113    int N;
    1214
    … …  
    1618    psLogSetFormat ("M");
    1719
    18     // these other options override the PSPHOT recipe options
     20    // these other options override the PSASTRO recipe options
    1921    psMetadata *options = psMetadataAlloc ();
    2022
    … …  
    2729
    2830    // load config data from default locations
    29     pmConfig *config = pmConfigRead(argc, argv);
     31    pmConfig *config = pmConfigRead (argc, argv);
    3032
    3133    // Storage for other command-line arguments
    … …  
    4345    }
    4446
    45     status = pmConfigFileSetsMD (config->arguments, "INPUT", "-file", "-list");
     47    status = pmConfigFileSetsMD (config->arguments, argc, argv, "INPUT", "-file", "-list");
    4648    if (!status) { usage ();}
    4749

  • trunk/psastro/src/psastroChipAstrom.c

    r6176 r7014  
    11# include "psastro.h"
    22
    3 // measure per-chip astrometry terms
    4 // this function requires a chip from a well-formed pmFPA
    5 // with PSASTRO.OBJECTS on the readouts
    6 bool psastroChipAstrom (pmChip *chip, psMetadata *config) {
     3bool psastroChipAstrom (pmConfig *config, psArray *refs) {
    74
    85    bool status;
    9     psArray *match;
    10     pmAstromStats stats;
     6    pmChip *chip = NULL;
     7    pmCell *cell = NULL;
     8    pmReadout *readout = NULL;
    119
    12     // cells->n > 1 is not well-defined
    13     if (chip->cells->n > 1) {
    14         psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nCells > 1");
    15         return false;
     10    // select the current recipe
     11    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO");
     12    if (!recipe) {
     13        psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n");
     14        exit(EXIT_FAILURE);
    1615    }
    17     pmCell *cell = chip->cells->data[0];
    1816
    19     // readouts->n > 1 is not well-defined
    20     if (cell->readouts->n > 1) {
    21         psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nReadouts > 1");
    22         return false;
     17    // select the input data sources
     18    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
     19    if (!input) {
     20        psErrorStackPrint(stderr, "Can't find input data!\n");
     21        exit(EXIT_FAILURE);
    2322    }
    24     pmReadout *readout = cell->readouts->data[0];
    25     psMetadata *header = pmReadoutGetHeader (readout);
    2623
    27     // load the corresponding reference data (DVO command)
    28     // XXX EAM : this needs to be improved, perhaps moved to a higher level
    29     psArray *refstars = psastroLoadReference (header, config);
     24    double plateScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLATE.SCALE");
     25    if (!status) plateScale = 1.0;
    3026
    31     // we require PASTRO.OBJECTS (pmAstromObj) on the analysis for this readout
    32     psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.OBJECTS");
     27    pmFPAview *view = pmFPAviewAlloc (0);
     28    pmFPA *fpa = input->fpa;
    3329
    34     // project the rawstars to the current best guess astrometry
    35     psastroProjectRawstars (rawstars, readout);
     30    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
     31        psTrace (__func__, 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
     32        if (!chip->process || !chip->file_exists) { continue; }
     33       
     34        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
     35            psTrace (__func__, 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
     36            if (!cell->process || !cell->file_exists) { continue; }
    3637
    37     // use the header & config info to project refstars onto the focal plane
    38     psastroProjectRefstars (refstars, readout);
     38            // process each of the readouts
     39            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
     40                if (! readout->data_exists) { continue; }
    3941
    40     // testWriteRaw ("ref.inp", refstars);
    41     // testWriteRaw ("raw.inp", rawstars);
     42                // select the raw objects for this readout
     43                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS");
     44                if (rawstars == NULL) { continue; }
    4245
    43     // fprintf (stderr, "rawstars:\n");
    44     // psastroDumpStars (rawstars);
    45     // fprintf (stderr, "refstars:\n");
    46     // psastroDumpStars (refstars);
     46                // the refstars is a subset within range of this chip
     47                psArray *refstars = psArrayAlloc (100);
    4748
    48     // find initial offset / rotation
    49     stats = pmAstromGridMatch (rawstars, refstars, config);
     49                // select the reference objects within range of this readout
     50                // project the reference objects to this chip
     51                for (int i = 0; i < refs->n; i++) {
     52                    pmAstromObj *ref = refs->data[i];
     53       
     54                    p_psProject (ref->TP, ref->sky, fpa->projection);
     55                    psPlaneDistortApply (ref->FP, fpa->fromTangentPlane, ref->TP, 0.0, 0.0);
     56                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
    5057
    51     // adjust the chip.toFPA terms only
    52     pmAstromGridApply (chip->toFPA, stats);
    53     chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
     58                    // XXX what is the X,Y range of the selected chip?
     59                    // XXX for the moment, force these to be good for Megacam
     60                    if (ref->chip->x < -500) continue;
     61                    if (ref->chip->x > 2500) continue;
     62                    if (ref->chip->y < -500) continue;
     63                    if (ref->chip->y > 5000) continue;
     64                   
     65                    psArrayAdd (refstars, 100, ref);
     66                }
     67                psastroOneChip (fpa, chip, refstars, rawstars, recipe);
     68                psFree (refstars);
    5469
    55     // use fit result to re-project rawstars
    56     psastroProjectRawstars (rawstars, readout);
    57     psastroProjectRefstars (refstars, readout);
    58 
    59     // testWriteRaw ("ref.dat", refstars);
    60     // testWriteRaw ("raw.dat", rawstars);
    61    
    62     // use small radius to match stars
    63     match = pmAstromRadiusMatch (rawstars, refstars, config);
    64 
    65     // improved fit for astrometric terms
    66     pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, config);
    67     chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
    68 
    69     // per-chip astrometry has been updated
    70     pmAstromWriteWCS (chip->toFPA, fpa->toSky, header);
    71 
     70                // read WCS data from the corresponding header
     71                pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
     72                pmAstromWriteWCS (chip->toFPA, fpa->toSky, hdu->header, plateScale);
     73            }
     74        }
     75    }
     76    psFree (view);
    7277    return true;
    7378}
     79
     80/* coordinate frame hierachy
     81   pixels (on a given readout)
     82   cell
     83   chip
     84   FP (focal plane)
     85   TP (tangent plane)
     86   sky (ra, dec)
     87*/

  • trunk/psastro/src/psastroConvert.c

    r6958 r7014  
    11# include "psastro.h"
     2// XXX leak free 2006.04.27
    23
     4// XXX pass/apply the WCS information?
     5bool psastroConvert (pmReadout *readout, psMetadata *recipe) {
     6
     7    psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
     8    if (sources == NULL)
     9        return false;
     10
     11    psArray *objects = pmSourceToAstromObj (sources);
     12
     13    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.OBJECTS", PS_DATA_ARRAY, "astrometry objects", objects);
     14    psFree (objects);
     15 
     16    return true;
     17}
     18
     19// XXX select a magnitude range?
    320psArray *pmSourceToAstromObj (psArray *sources) {
    421
    … …  
    623   
    724    for (int i = 0; i < sources->n; i++) {
    8         psSource *source = sources->data[i];
     25        pmSource *source = sources->data[i];
    926       
    1027        // only accept the PSF sources?
    11         if (source->type != PM_SOURCE_STAR) continue;
     28        if (source->type != PM_SOURCE_TYPE_STAR) continue;
    1229
    1330        pmModel *model = source->modelPSF;
    … …  
    2542        // XXX do we have the information giving the readout and cell offset?
    2643        // for the moment, assume chip == cell == readout
    27         obj->cell = obj->pix;
    28         obj->chip = obj->cell;
    29        
    30         // can we make an initial guess of the astrometry to set FP, TP, and sky coords?
    31        
    32        
     44        *obj->cell = *obj->pix;
     45        *obj->chip = *obj->cell;
     46
     47        psArrayAdd (objects, 100, obj);
     48        psFree (obj);
     49    }
     50    return objects;
    3351}

  • trunk/psastro/src/psastroDataSave.c

    r6792 r7014  
    11# include "psastro.h"
     2// XXX leak free 2006.04.27
    23
    3 // load the data from the files in this loop.
    4 // we write out the result in a second loop
    5 // at the end of this function, the complete stellar data is loaded
    6 // into the correct fpa structure locations (readout.analysis:PSPHOT.SOURCES)
     4// this loop saves the photometry/astrometry data files
    75bool psastroDataSave (pmConfig *config) {
    86
    9   // define the output files:
    10     pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSASTRO.INPUT");
    11     if (!status) {
    12         psErrorStackPrint(stderr, "Can't find input data!\n");
     7    pmChip *chip;
     8    pmCell *cell;
     9    pmReadout *readout;
     10
     11    // select the current recipe
     12    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO");
     13    if (!recipe) {
     14        psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n");
    1315        exit(EXIT_FAILURE);
    1416    }
    1517
    16     // this is the output target
    17     pmFPAfileDefine (config->files, format, input, "PSASTRO.OUTPUT");
     18    // select the output data sources
     19    pmFPAfile *output = psMetadataLookupPtr (NULL, config->files, "PSASTRO.OUTPUT");
     20    if (!output) {
     21        psErrorStackPrint(stderr, "Can't find output data!\n");
     22        exit(EXIT_FAILURE);
     23    }
    1824
    19     char *output = psMetadataLookupPtr(&status, config->arguments, "OUTPUT");
    20     pmFPAfileAddOutputName (config->files, "OUTPUT", output, PM_FPA_MODE_WRITE);
     25    // de-activate all files except PSASTRO.OUTPUT
     26    pmFPAfileActivate (config->files, false, NULL);
     27    pmFPAfileActivate (config->files, true, "PSASTRO.OUTPUT");
    2128
    2229    pmFPAview *view = pmFPAviewAlloc (0);
    2330
    24     // files associated with the science image
     31    // open/load files as needed
    2532    pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE);
    2633
    27     while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    28         psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
     34    while ((chip = pmFPAviewNextChip (view, output->fpa, 1)) != NULL) {
     35        psTrace (__func__, 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    2936        if (!chip->process || !chip->file_exists) { continue; }
    3037        pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE);
    3138
    32         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
    33             psLogMsg ("psphot", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
     39        while ((cell = pmFPAviewNextCell (view, output->fpa, 1)) != NULL) {
     40            psTrace (__func__, 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    3441            if (!cell->process || !cell->file_exists) { continue; }
    3542            pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE);
    3643
    3744            // process each of the readouts
    38             while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
     45            while ((readout = pmFPAviewNextReadout (view, output->fpa, 1)) != NULL) {
    3946                pmFPAfileIOChecks (config->files, view, PM_FPA_BEFORE);
    4047                if (! readout->data_exists) { continue; }
    … …  
    4754    }
    4855    pmFPAfileIOChecks (config->files, view, PM_FPA_AFTER);
     56
     57    psFree (view);
    4958    return true;
    5059}

  • trunk/psastro/src/psastroParseCamera.c

    r6911 r7014  
    1 # include "psphot.h"
     1# include "psastro.h"
     2// XXX leak free 2006.04.27
    23
    3 bool *psastroParseCamera (pmConfig *config) {
     4bool psastroParseCamera (pmConfig *config) {
    45
    5     ppFile *input = ppFileAlloc ();
     6    bool status = false;
    67
    78    // the input image(s) are required arguments; they define the camera
    8     input = pmFPAfileAddSource (config, "INPUT", "PSASTRO.INPUT", true);
     9    pmFPAfile *input = pmFPAfileFromArgs (&status, config, "PSASTRO.INPUT", "INPUT");
    910    if (!status) { psAbort (__func__, "missing INPUT entry"); }
    1011
    … …  
    1516
    1617    // set default recipe values here
    17     psMetadataAddStr (recipe, PS_LIST_TAIL, "PHOTCODE", PS_META_NO_REPLACE, "default fitting mode", "NONE");
    18     psMetadataAddStr (recipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_NO_REPLACE, "default fitting mode", "NONE");
     18    psMetadataAddStr (recipe, PS_LIST_TAIL, "PHOTCODE",    PS_META_NO_REPLACE, "", "NONE");
     19    psMetadataAddStr (recipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_NO_REPLACE, "", "NONE");
    1920
    2021    // these calls bind the I/O handle to the specified fpa
    … …  
    3839        }
    3940    }
     41    psFree (chips);
    4042
    4143    psTrace(__func__, 1, "Done with psastroParseCamera...\n");
    4244    return true;
    4345}
     46
     47
     48// useful for debugging
     49# if 0
     50for (int i = 0; i < input->fpa->chips->n; i++) {
     51    pmChip *chip = input->fpa->chips->data[i];
     52    fprintf (stderr, "chip %2d: %x %x\n", i, chip->file_exists, chip->process);
     53
     54    for (int j = 0; j < chip->cells->n; j++) {
     55        pmCell *cell = chip->cells->data[j];
     56        fprintf (stderr, "cell %2d: %x %x\n", j, cell->file_exists, cell->process);
     57
     58    }
     59}
     60# endif

  • trunk/psastro/src/psastroUtils.c

    r6176 r7014  
    11# include "psastro.h"
    22# define RENORM 0
     3
     4bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars) {
     5
     6    // XXX this region needs to be defined more intelligently
     7    psRegion region = {0, 0, 2000, 4500};
     8    psFree (chip->fromFPA);
     9    chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, region, 20);
     10
     11    for (int i = 0; i < rawstars->n; i++) {
     12        pmAstromObj *raw = rawstars->data[i];
     13       
     14        psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
     15        psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0);
     16        p_psDeproject (raw->sky, raw->TP, fpa->projection);
     17    }
     18
     19    for (int i = 0; i < refstars->n; i++) {
     20        pmAstromObj *ref = refstars->data[i];
     21        psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
     22    }
     23
     24    return true;
     25}
     26
     27# if 0
    328
    429bool psastroSelectBrightStars (pmFPA *fpa, psMetadata *config) {
    … …  
    143168    return true;
    144169}
     170
     171psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input) {
     172
     173    psPolynomial2D *output = psPolynomial2DAlloc (input->nX, input->nY, input->type);
     174
     175    for (int i = 0; i < input->nX; i++) {
     176        for (int j = 0; j < input->nY; j++) {
     177            output->mask[i][j]     = input->mask[i][j];
     178            output->coeff[i][j]    = input->coeff[i][j];
     179            output->coeffErr[i][j] = input->coeffErr[i][j];
     180        }
     181    }
     182    return (output);
     183}
     184
     185psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input) {
     186
     187    psPolynomial4D *output = psPolynomial4DAlloc (input->nX, input->nY, input->nZ, input->nT, input->type);
     188
     189    for (int i = 0; i < input->nX; i++) {
     190        for (int j = 0; j < input->nY; j++) {
     191            for (int k = 0; k < input->nZ; k++) {
     192                for (int m = 0; m < input->nT; m++) {
     193                    output->mask[i][j][k][m]     = input->mask[i][j][k][m];
     194                    output->coeff[i][j][k][m]    = input->coeff[i][j][k][m];
     195                    output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m];
     196                }
     197            }
     198        }
     199    }
     200    return (output);
     201}
     202
     203psPlaneDistort *psPlaneDistortCopy (psPlaneDistort *input) {
     204
     205    psPlaneDistort *output = psPlaneDistortAlloc (input->x->nX, input->x->nY, input->x->nZ, input->x->nT);
     206
     207    for (int i = 0; i < input->x->nX; i++) {
     208        for (int j = 0; j < input->x->nY; j++) {
     209            for (int k = 0; k < input->x->nZ; k++) {
     210                for (int m = 0; m < input->x->nT; m++) {
     211                    // x-terms
     212                    output->x->mask[i][j][k][m]     = input->x->mask[i][j][k][m];
     213                    output->x->coeff[i][j][k][m]    = input->x->coeff[i][j][k][m];
     214                    output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m];
     215                    // y-terms
     216                    output->y->mask[i][j][k][m]     = input->y->mask[i][j][k][m];
     217                    output->y->coeff[i][j][k][m]    = input->y->coeff[i][j][k][m];
     218                    output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m];
     219                }
     220            }
     221        }
     222    }
     223    return (output);
     224}
     225
     226psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input) {
     227
     228    psPlaneTransform *output = psPlaneTransformAlloc (input->x->nX, input->x->nY);
     229
     230    for (int i = 0; i < input->x->nX; i++) {
     231        for (int j = 0; j < input->x->nY; j++) {
     232            // x-terms
     233            output->x->mask[i][j]     = input->x->mask[i][j];
     234            output->x->coeff[i][j]    = input->x->coeff[i][j];
     235            output->x->coeffErr[i][j] = input->x->coeffErr[i][j];
     236            // y-terms
     237            output->y->mask[i][j]     = input->y->mask[i][j];
     238            output->y->coeff[i][j]    = input->y->coeff[i][j];
     239            output->y->coeffErr[i][j] = input->y->coeffErr[i][j];
     240        }
     241    }
     242    return (output);
     243}
     244
     245psProjection *psProjectionCopy (psProjection *input) {
     246
     247    psProjection *output = psProjectionAlloc (input->R, input->D, input->Xs, input->Ys, input->type);
     248    return (output);
     249}
     250
     251// very crude distortion inversion: assumes 0 order in z and t, linear in x and y:
     252psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort) {
     253    PS_ASSERT_PTR_NON_NULL(distort, 0);
     254    PS_ASSERT_PTR_NON_NULL(distort->x, 0);
     255    PS_ASSERT_PTR_NON_NULL(distort->y, 0);
     256
     257    psPlaneDistort *out = psPlaneDistortAlloc(1, 1, 0, 0);
     258
     259    /* simple matrix inversion code */
     260
     261    psF64 r11 = distort->x->coeff[1][0][0][0];
     262    psF64 r12 = distort->x->coeff[0][1][0][0];
     263    psF64 r21 = distort->y->coeff[1][0][0][0];
     264    psF64 r22 = distort->y->coeff[0][1][0][0];
     265    psF64 xo  = distort->x->coeff[0][0][0][0];
     266    psF64 yo  = distort->y->coeff[0][0][0][0];
     267
     268    psF64 invDet = 1.0 / (r11 * r22 - r12 * r21); // Inverse of the determinant
     269
     270    out->x->coeff[1][0][0][0] = +invDet * r22;
     271    out->x->coeff[0][1][0][0] = -invDet * r12;
     272    out->y->coeff[1][0][0][0] = -invDet * r21;
     273    out->y->coeff[0][1][0][0] = +invDet * r11;
     274
     275    out->x->coeff[0][0][0][0] = - invDet * (r22 * xo - r12 * yo);
     276    out->y->coeff[0][0][0][0] = - invDet * (r11 * yo - r21 * xo);
     277
     278    return(out);
     279}
     280
     281// returns the rotation term, forcing positive parity
     282double psPlaneTransformGetRotation (psPlaneTransform *map) {
     283
     284    if (map->x->nX < 1) return 0;
     285    if (map->x->nY < 1) return 0;
     286
     287    if (map->y->nX < 1) return 0;
     288    if (map->y->nY < 1) return 0;
     289   
     290    double pc1_1 = map->x->coeff[1][0];
     291    double pc1_2 = map->x->coeff[0][1];
     292    double pc2_1 = map->y->coeff[1][0];
     293    double pc2_2 = map->y->coeff[0][1];
     294
     295    double px = SIGN (pc1_1);
     296    double py = SIGN (pc2_2);
     297
     298    // both x and y terms imply an angle. take the average
     299    double t1 = -atan2 (px*pc1_2, px*pc1_1);
     300    double t2 = +atan2 (py*pc2_1, py*pc2_2);
     301   
     302    // careful near -pi,+pi boundary...
     303    if (t1 - t2 > M_PI/2) t2 += 2*M_PI;
     304    if (t2 - t1 > M_PI/2) t1 += 2*M_PI;
     305
     306    double theta = 0.5*(t1 + t2);
     307    while (theta < M_PI) theta += 2*M_PI;
     308    while (theta > M_PI) theta -= 2*M_PI;
     309   
     310    return (theta);
     311}
     312
     313# endif

  • trunk/psastro/src/psastroWCS.c

    r6176 r7014  
    22
    33// interpret header WCS (only handles traditional WCS for the moment)
    4 bool pmAstromReadWCS (psPlaneTransform **toFPAOut, psProjection **toSkyOut, psMetadata *header) {
     4// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
     5bool pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale) {
    56
    67    psPlaneTransform *toFPA;
     8    psPlaneDistort *toTPA;
    79    psProjection *toSky;
    810    psProjectionType type;
    … …  
    1921
    2022    // determine projection type
     23    // XXX add the Elixir DIS / WRP two-layer projection here
    2124    type = PS_PROJ_NTYPE;
    2225    if (!strcmp (&ctype[4], "-SIN")) type = PS_PROJ_SIN;
    … …  
    7275       
    7376        // renormalize to cdelt1, cdelt2, etc
    74         float scale = hypot (pc1_1, pc1_2);
     77        double scale = hypot (pc1_1, pc1_2);
    7578        cdelt1 = cdelt2 = scale;
    7679        pc1_1 /= scale;
    … …  
    8992    //           scale = scale(i)/scale(0) (i == chip #)
    9093    //           project crval1(0),crval2(0 using projection
    91 
    92     // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
    93     // XXX EAM : I've fixed this in pslib eam_rel8_b2
    9494    toFPA = psPlaneTransformAlloc (1, 1);
    9595   
    96     toFPA->x->coeff[0][0] = crpix1;
    97     toFPA->x->coeff[1][0] = pc1_1;
    98     toFPA->x->coeff[0][1] = pc1_2;
     96    double cdelt = hypot (cdelt1, cdelt2) / plateScale;  // degrees / micron (eg, in fact, whatever unit user chooses for focal plane)
     97    cdelt1 /= cdelt;
     98    cdelt2 /= cdelt;
     99
     100    toFPA->x->coeff[0][0] = -(pc1_1*cdelt1*crpix1 + pc1_2*cdelt2*crpix2);
     101    toFPA->x->coeff[1][0] = +(pc1_1*cdelt1);
     102    toFPA->x->coeff[0][1] = +(pc1_2*cdelt2);
    99103    toFPA->x->mask[1][1]  = 1;
    100104
    101     toFPA->y->coeff[0][0] = crpix2;
    102     toFPA->y->coeff[1][0] = pc2_1;
    103     toFPA->y->coeff[0][1] = pc2_2;
     105    toFPA->y->coeff[0][0] = -(pc2_1*cdelt1*crpix1 + pc2_2*cdelt2*crpix2);
     106    toFPA->y->coeff[1][0] = +(pc2_1*cdelt1);
     107    toFPA->y->coeff[0][1] = +(pc2_2*cdelt2);
    104108    toFPA->y->mask[1][1]  = 1;
    105 
    106     // center of projection is (0,0) coordinate of TPA
    107     toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
    108 
    109     *toFPAOut = toFPA;
    110     *toSkyOut = toSky;
    111 
     109    *toFPAout = toFPA;
     110
     111    if (*toSkyOut != NULL) {
     112        if (*toTPAout == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
     113
     114        psSphere sky;
     115        psPlane tpa;
     116
     117        sky.r = crval1*RAD_DEG;
     118        sky.d = crval2*RAD_DEG;
     119        p_psProject (&tpa, &sky, *toSkyOut);
     120       
     121        // XXX for the moment, assume toTPA is the identity transformation
     122        toFPA->x->coeff[0][0] = tpa.x;
     123        toFPA->y->coeff[0][0] = tpa.y;
     124    } else {
     125        // XXX for now, use the identity for TPA <--> FPA
     126        toTPA = psPlaneDistortIdentity ();
     127
     128        // center of projection is (0,0) coordinate of TPA
     129        toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, RAD_DEG*cdelt, RAD_DEG*cdelt, type);
     130
     131        *toTPAout = toTPA;
     132        *toSkyOut = toSky;
     133    }
    112134    return true;
    113135}
    … …  
    115137
    116138// convert toFPA / toSky components to traditional WCS
    117 bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header) {
     139// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
     140bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale) {
    118141
    119142    switch (toSky->type) {
    120143      case PS_PROJ_SIN:
    121         psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---SIN");
    122         psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--SIN");
     144        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---SIN");
     145        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--SIN");
    123146        break;
    124147      case PS_PROJ_TAN:
    125         psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---TAN");
    126         psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--TAN");
     148        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---TAN");
     149        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--TAN");
    127150        break;
    128151      case PS_PROJ_AIT:
    129         psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---AIT");
    130         psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--AIT");
     152        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---AIT");
     153        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--AIT");
    131154        break;
    132155      case PS_PROJ_PAR:
    133         psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---PAR");
    134         psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--PAR");
     156        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---PAR");
     157        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--PAR");
    135158        break;
    136159      default:
    … …  
    139162    }
    140163
    141     psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->R*DEG_RAD);
    142     psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->D*DEG_RAD);
    143     psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
    144     psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
    145     psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Xs*DEG_RAD);
    146     psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Ys*DEG_RAD);
    147 
    148     psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[1][0]);
    149     psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][1]);
    150     psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[1][0]);
    151     psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][1]);
     164    // XXX not really right: needs to deal with non-identity toTP coeffs
     165    // XXX actually, totally wrong.  fix the conversions
     166    // XXX need to handle the plateScale
     167    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",   PS_META_REPLACE, "", toSky->R*DEG_RAD);
     168    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",   PS_META_REPLACE, "", toSky->D*DEG_RAD);
     169    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",   PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
     170    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",   PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
     171    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",   PS_META_REPLACE, "", toSky->Xs*DEG_RAD);
     172    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",   PS_META_REPLACE, "", toSky->Ys*DEG_RAD);
     173
     174    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]);
     175    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]);
     176    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]);
     177    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]);
    152178
    153179    // alternative representations use
    … …  
    157183    return true;
    158184}
     185
     186psPlaneDistort *psPlaneDistortIdentity () {
     187
     188    psPlaneDistort *distort;
     189
     190    distort = psPlaneDistortAlloc (1, 1, 0, 0);
     191    distort->x->coeff[0][0][0][0] = 0;
     192    distort->x->coeff[1][0][0][0] = 1;
     193    distort->x->coeff[0][1][0][0] = 0;
     194    distort->x->mask [1][1][0][0] = 1;
     195
     196    distort->y->coeff[0][0][0][0] = 0;
     197    distort->y->coeff[1][0][0][0] = 0;
     198    distort->y->coeff[0][1][0][0] = 1;
     199    distort->y->mask [1][1][0][0] = 1;
     200
     201    return distort;
     202}
Note: See TracChangeset for help on using the changeset viewer.