IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 23335


Ignore:
Timestamp:
Mar 15, 2009, 4:56:07 PM (17 years ago)
Author:
eugene
Message:

add psastroExtract (first working version)

Location:
branches/eam_branches/eam_branch_20090312/psastro/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtract.c

    r23322 r23335  
    4141    }
    4242
    43     // load the reference stars overlapping the data stars
    44     psArray *refs = psastroLoadRefstars(config);
    45     if (!refs) {
    46         psError (PSASTRO_ERR_UNKNOWN, false, "failed to load reference data\n");
    47         return false;
    48     }
    49     if (refs->n == 0) {
    50         psError(PSASTRO_ERR_REFSTARS, true, "no reference stars found");
    51         psFree(refs);
    52         return false;
    53     }
    54 
    55     if (!psastroChooseRefstars (config, refs)) {
    56         psError (PSASTRO_ERR_UNKNOWN, false, "failed to select reference data for chips\n");
    57         psFree(refs);
    58         return false;
    59     }
    60 
    6143    // run the full astrometry analysis (chip and/or mosaic)
    6244    if (!psastroExtractAnalysis (config)) {
    … …  
    6547    }
    6648   
    67     // run the full astrometry analysis (chip and/or mosaic)
    68     if (!psastroExtractAdjust (config)) {
    69         psErrorStackPrint(stderr, "failure in psastro model adjust\n");
    70         exit (1);
    71     }
    72    
    73     // save the model
    74     if (!psastroExtractDataSave (config)) {
    75         psErrorStackPrint(stderr, "error saving output data\n");
    76         exit (1);
    77     }
     49    // XXX save the results; is this needed?
     50    // if (!psastroExtractDataSave (config)) {
     51    // psErrorStackPrint(stderr, "error saving output data\n");
     52    // exit (1);
     53    // }
    7854
    7955    psLogMsg ("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete"));
    8056
    81     // psastroCleanup (config);
     57    psastroCleanup (config);
    8258    exit (EXIT_SUCCESS);
    8359}

  • branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractAnalysis.c

    r23322 r23335  
    2020
    2121# define TEST_OUTPUT 1
    22 static pmFits *outStars = NULL;
     22
     23psImage *psastroExtractStar (psImage *input, float x, float y, float dX, float dY);
    2324
    2425/**
    … …  
    3233    pmReadout *readout = NULL;
    3334    float zeropt, exptime;
     35    char extname[81];
    3436
    3537    // select the current recipe
    … …  
    4042    }
    4143
     44    // force this recipe value to be false (would cause an error in psastroChooseRefstars -- no rawstars!)
     45    psMetadataAddBool (recipe, PS_LIST_TAIL, "PSASTRO.MATCH.LUMFUNC", PS_META_REPLACE, "do not match luminosity functions", false);
     46
    4247    psLogMsg ("psastro", PS_LOG_INFO, "extracting bright-star subrasters");
     48
     49    // load the reference stars overlapping the data stars
     50    psArray *refs = psastroLoadRefstars(config, "PSASTRO.EXTRACT.ASTROM");
     51    if (!refs) {
     52        psError (PSASTRO_ERR_UNKNOWN, false, "failed to load reference data\n");
     53        return false;
     54    }
     55    if (refs->n == 0) {
     56        psError(PSASTRO_ERR_REFSTARS, true, "no reference stars found");
     57        psFree(refs);
     58        return false;
     59    }
     60
     61    if (!psastroChooseRefstars (config, refs, "PSASTRO.EXTRACT.ASTROM")) {
     62        psError (PSASTRO_ERR_UNKNOWN, false, "failed to select reference data for chips\n");
     63        psFree(refs);
     64        return false;
     65    }
     66    psFree(refs);
    4367
    4468    double EXTRACT_MAX_MAG = psMetadataLookupF32 (&status, recipe, "EXTRACT_MAX_MAG");
    … …  
    4670    // we have two input pmFPAfiles: PSASTRO.EXTRACT.INPUT and PSASTRO.EXTRACT.ASTROM.  both are in chip-mosaic or fpa format
    4771    // we have already loaded the headers from PSASTRO.EXTRACT.ASTROM and determined the astrometry terms
     72
     73    // select the input astrometry data (also carries the refstars)
     74    pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.EXTRACT.ASTROM");
     75    if (!astrom) {
     76        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     77        return false;
     78    }
     79    pmFPA *fpa = astrom->fpa;
    4880
    4981    // select the input data sources
    … …  
    5385        return false;
    5486    }
    55     pmFPA *fpa = input->fpa;
    56 
    57     // select the input data sources
    58     pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.EXTRACT.ASTROM");
    59     if (!input) {
    60         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    61         return false;
    62     }
    6387
    6488    // only load data from INPUT
    … …  
    77101    EXTRACT_MAX_MAG += MagOffset;
    78102
     103    int nExt = 0;
     104
    79105    // open the output file handle: we are just saving a series of extensions to this file
     106    char *filename = psMetadataLookupStr (&status, config->arguments, "OUTPUT");
     107    if (!filename) {
     108      psLogMsg ("psastro", PS_LOG_INFO, "output filename not supplieda");
     109      return false;
     110    }
     111    psFits *outStars = psFitsOpen (filename, "w");
    80112    if (!outStars) {
    81       char *filename = psMetadataLookupStr (&stats, config->arguments, "OUTPUT");
    82       if (!filename) {
    83         psLogMsg ("psastro", PS_LOG_INFO, "output filename not supplieda");
    84         return false;
    85       }
    86       outStars = psFitsOpen (filename, "w");
    87       if (!outStars) {
    88         psError(PS_ERR_IO, false, "error opening file %s\n", filename);
    89         return false;
    90       }
    91       // call: psFitsWriteImage (fits, NULL, subraster, 0, "extname.NN");
    92     }
     113      psError(PS_ERR_IO, false, "error opening file %s\n", filename);
     114      return false;
     115    }
     116
     117    pmFPAview *view = pmFPAviewAlloc (0);
    93118
    94119    // open/load files as needed
    … …  
    109134          char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME");
    110135          psStringAppend (&filename, "bright.%s.dat", chipname);
    111           FILE *f = fopen (filename, "w");
     136          f = fopen (filename, "w");
    112137          if (!f) {
    113138            psWarning ("cannot create test output file %s\n", filename);
    … …  
    132157                for (int i = 0; i < refstars->n; i++) {
    133158                    pmAstromObj *ref = refstars->data[i];
    134                     if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue;
     159                    if (ref->Mag > EXTRACT_MAX_MAG) continue;
    135160
    136161                    if (TEST_OUTPUT) {
    … …  
    138163                    }
    139164
    140                     // XXX for now, assume cell binning is 1x1 relative to chip
    141                     psImage *subraster = psastroExtractStar (input->image, ref->chip->x, ref->chip->y, DX, DY);
     165                    pmReadout *inReadout = pmFPAviewThisReadout (view, input->fpa);
     166                    psImage *subraster = psastroExtractStar (inReadout->image, ref->chip->x, ref->chip->y, 100.0, 100.0);
     167                    if (!subraster) continue;
     168                   
     169                    snprintf (extname, 80, "extname.%05d", nExt);
     170                    psFitsWriteImage (outStars, NULL, subraster, 0, extname);
     171                    nExt ++;
     172                   
     173                    psFree (subraster);
    142174                }
    143                 if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE;
     175                if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
    144176            }
    145             if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE;
     177            if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
    146178        }
    147179        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
    148180
     181        fclose (f);
    149182    }
    150183    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
    151184
    152     // need to close the output fits file on some basis..
    153     fclose (f);
    154 
     185    psFitsClose (outStars);
    155186    psFree (view);
    156187    return true;
    157188}
    158189
     190psImage *psastroExtractStar (psImage *input, float x, float y, float dX, float dY) {
     191
     192  // skip if no pixels are on the image
     193  // skip if center is not on the image
     194  // if bounds fall off image, paste into a full-size image.
     195
     196  if (x < 0) return NULL;
     197  if (y < 0) return NULL;
     198  if (x >= input->numCols) return NULL;
     199  if (y >= input->numRows) return NULL;
     200
     201  psRegion fullRegion = psRegionSet (x - dX, x + dX, y - dY, y + dY);
     202  psRegion realRegion = psRegionForImage (input, fullRegion);
     203  psImage *subset = psImageSubset (input, realRegion);
     204
     205  // if ((subset->numCols < dX) || (subset->numRows < dY)) {
     206  //   psImage *fullsize = psImageAlloc (dX, dY, PS_TYPE_F32);
     207  //   // fill in the fullsize image with the values from the subset
     208  //   // I must already have done this elsewhere...
     209  // }
     210  return subset;
     211}

  • branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractArguments.c

    r23322 r23335  
    3737    // psMetadata *options = pmConfigRecipeOptions (config, PSASTRO_RECIPE);
    3838
     39    // define the image pixel data
    3940    if (!pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list")) {
    4041        psError(PSASTRO_ERR_ARGUMENTS, true, "Missing -file (input) or -list (input)");
    … …  
    4344    }
    4445
    45     // XXX require this?
    46     pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM",   "-astrom", "-astromlist");
     46    // define the astrometry data (XXX make this optional and use the image header?)
     47    if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM",   "-astrom", "-astromlist")) {
     48        psError(PSASTRO_ERR_ARGUMENTS, true, "Missing -astrom (smf) or -astromlist (smf)");
     49        psErrorStackPrint(stderr, "%s", usage);
     50        exit (1);
     51    }
    4752
    4853    // define the output filename

  • branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractDataLoad.c

    r23322 r23335  
    2929    pmChip *chip;
    3030
     31    bool newFPA = true;
     32    double RAmin  = +FLT_MAX;
     33    double RAmax  = -FLT_MAX;
     34    double DECmin = +FLT_MAX;
     35    double DECmax = -FLT_MAX;
     36
     37    double RAminSky = NAN;
     38    double RAmaxSky = NAN;
     39
    3140    psTimerStart ("psastro");
    3241
    … …  
    4857    pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.EXTRACT.ASTROM");
    4958    if (!astrom) psAbort ("PSASTRO.EXTRACT.ASTROM not listed in config->files");
     59    pmFPA *fpa = astrom->fpa;
    5060
    5161    pmFPAfileActivate(config->files, false, NULL);
    … …  
    7686    }
    7787
     88    // apply the header astrometry to the astrometry structures
    7889    while ((chip = pmFPAviewNextChip (view, astrom->fpa, 1)) != NULL) {
    7990      psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    8091      if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
     92
     93      if (newFPA) {
     94        newFPA = false;
     95        while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
     96        while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
     97        RAminSky = fpa->toSky->R - M_PI;
     98        RAmaxSky = fpa->toSky->R + M_PI;
     99      }
    81100
    82101      // read WCS data from the corresponding header
    … …  
    96115        }
    97116      }
     117     
     118      // determine RA,DEC of 4 corners, use to find RA_MIN,MAX, DEC_MIN,MAX
     119      psRegion *region = pmChipPixels (chip);
     120      psPlane ptCH[4], ptFP, ptTP;
     121      psSphere ptSky;
     122      ptCH[0].x = region->x0;
     123      ptCH[0].y = region->y0;
     124      ptCH[1].x = region->x1;
     125      ptCH[1].y = region->y0;
     126      ptCH[2].x = region->x1;
     127      ptCH[2].y = region->y1;
     128      ptCH[3].x = region->x0;
     129      ptCH[3].y = region->y1;
     130      psFree (region);
     131     
     132      for (int i = 0; i < 4; i++) {
     133        psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH[i]);
     134        psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
     135        psDeproject (&ptSky, &ptTP, fpa->toSky);
     136
     137        // rationalize ra to sky range centered on boresite
     138        while (ptSky.r < RAminSky) ptSky.r += 2.0*M_PI;
     139        while (ptSky.r > RAmaxSky) ptSky.r -= 2.0*M_PI;
     140
     141        RAmin = PS_MIN (ptSky.r, RAmin);
     142        RAmax = PS_MAX (ptSky.r, RAmax);
     143       
     144        DECmin = PS_MIN (ptSky.d, DECmin);
     145        DECmax = PS_MAX (ptSky.d, DECmax);
     146
     147        psLogMsg ("psastro", 2, "chip %d (corner %d) : %f %f  -> %f %f -> %f %f -> %f %f\n",
     148                  view->chip, i, ptCH[i].x, ptCH[i].y, ptFP.x, ptFP.y, ptTP.x, ptTP.y,
     149                  DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
     150      }
    98151    }
    99152    psLogMsg ("psastro", 3, "convert wcs terms to internal format : %f sec\n", psTimerMark ("psastro"));
     153
     154    psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n",
     155              DEG_RAD*RAmin, DEG_RAD*DECmin,
     156              DEG_RAD*RAmax, DEG_RAD*DECmax);
     157
     158    psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN",  PS_META_REPLACE, "", RAmin);
     159    psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MAX",  PS_META_REPLACE, "", RAmax);
     160    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DEC_MIN", PS_META_REPLACE, "", DECmin);
     161    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DEC_MAX", PS_META_REPLACE, "", DECmax);
    100162
    101163    psFree (view);

  • branches/eam_branches/eam_branch_20090312/psastro/src/psastroExtractParseCamera.c

    r23322 r23335  
    1616
    1717    bool status = false;
    18     pmFPAfile *input = NULL;
    1918
    2019    // the input image(s) are required arguments; they define the camera
    … …  
    2928    status = false;
    3029    pmFPAfile *astrom = pmFPAfileDefineFromArgs(&status, config, "PSASTRO.EXTRACT.ASTROM", "ASTROM");
    31     if (!status) {
     30    if (!astrom || !status) {
    3231        psError (PS_ERR_UNKNOWN, false, "failed to load astrometry definition");
    3332        return NULL;
    34     }
    35     if (astrom) {
    36         psLogMsg ("psastro", 3, "using supplied astrometry\n");
    37     } else {
    38         psLogMsg ("psastro", 3, "using header astrometry\n");
    3933    }
    4034
    … …  
    5448    psArray *chips = psStringSplitArray (chipLine, ",", false);
    5549    if (chips->n > 0) {
    56         pmFPASelectChip (output->fpa, -1, true); // deselect all chips
     50        pmFPASelectChip (input->fpa, -1, true); // deselect all chips
     51        pmFPASelectChip (astrom->fpa, -1, true); // deselect all chips
    5752        for (int i = 0; i < chips->n; i++) {
    5853            int chipNum = atoi(chips->data[i]);
    59             if (! pmFPASelectChip(output->fpa, chipNum, false)) {
     54            if (! pmFPASelectChip(input->fpa, chipNum, false)) {
     55                psError(PSASTRO_ERR_CONFIG, true, "Chip number %d doesn't exist in camera.\n", chipNum);
     56                return false;
     57            }
     58            if (! pmFPASelectChip(astrom->fpa, chipNum, false)) {
    6059                psError(PSASTRO_ERR_CONFIG, true, "Chip number %d doesn't exist in camera.\n", chipNum);
    6160                return false;
Note: See TracChangeset for help on using the changeset viewer.