IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15359


Ignore:
Timestamp:
Oct 23, 2007, 10:54:53 AM (19 years ago)
Author:
eugene
Message:

adding pieces for astrometry table I/O, extended source analysis

Location:
branches/eam_branch_20071023
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20071023/psModules/src/camera/pmFPAfileIO.c

    r15333 r15359  
    192192        status = pmPSFmodelReadForView (view, file, config);
    193193        break;
     194      case PM_FPA_FILE_ASTROM:
     195        status = pmAstromReadForView (view, file, config);
     196        break;
    194197      case PM_FPA_FILE_JPEG:
    195198      case PM_FPA_FILE_KAPA:
    … …  
    273276    case PM_FPA_FILE_CMF:
    274277    case PM_FPA_FILE_PSF:
     278    case PM_FPA_FILE_ASTROM:
    275279    case PM_FPA_FILE_JPEG:
    276280    case PM_FPA_FILE_KAPA:
    … …  
    336340    // have their data stored on the readout->analysis metadata structure of another
    337341    // (existing) fpa
     342    if (file->type == PM_FPA_FILE_CMF) {
     343      if (!pmFPAviewCheckDataStatusForSources (view, file)) {
     344        psTrace("psModules.camera", 6, "skip write for %s, no data for this entry", file->name);
     345        return true;
     346      }
     347    }
    338348    if (file->type == PM_FPA_FILE_PSF) {
    339349      if (!pmPSFmodelCheckDataStatusForView (view, file)) {
    … …  
    342352      }
    343353    }
    344     if (file->type == PM_FPA_FILE_CMF) {
    345       if (!pmFPAviewCheckDataStatusForSources (view, file)) {
     354    if (file->type == PM_FPA_FILE_ASTROM) {
     355      if (!pmAstromCheckDataStatusForView (view, file)) {
    346356        psTrace("psModules.camera", 6, "skip write for %s, no data for this entry", file->name);
    347357        return true;
    … …  
    449459      case PM_FPA_FILE_PSF:
    450460        status = pmPSFmodelWriteForView (view, file, config);
     461        break;
     462
     463      case PM_FPA_FILE_ASTROM:
     464        status = pmAstromWriteForView (view, file, config);
    451465        break;
    452466
    … …  
    508522      case PM_FPA_FILE_CMF:
    509523      case PM_FPA_FILE_PSF:
     524      case PM_FPA_FILE_ASTROM:
    510525        psTrace ("psModules.camera", 5, "closing %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout);
    511526        status = psFitsClose (file->fits);
    … …  
    574589      case PM_FPA_FILE_CMF:
    575590      case PM_FPA_FILE_PSF:
     591      case PM_FPA_FILE_ASTROM:
    576592        psTrace ("psModules.camera", 5, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name);
    577593        return true;
    … …  
    712728      case PM_FPA_FILE_CMF:
    713729      case PM_FPA_FILE_PSF:
     730      case PM_FPA_FILE_ASTROM:
    714731        psTrace ("psModules.camera", 5, "opening %s (%s) (%d:%d:%d)\n",
    715732                 file->filename, file->name, view->chip, view->cell, view->readout);
    … …  
    861878        status = pmPSFmodelWritePHU (view, file, config);
    862879        break;
     880      case PM_FPA_FILE_ASTROM:
     881        status = pmAstromWritePHU (view, file, config);
     882        break;
    863883      case PM_FPA_FILE_SX:
    864884      case PM_FPA_FILE_RAW:

  • branches/eam_branch_20071023/psModules/src/objects/pmModel.c

    r14938 r15359  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2007-09-21 00:04:07 $
     8 *  @version $Revision: 1.15.6.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2007-10-23 20:54:53 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    217217        Ro = psImageInterpolateOptionsAlloc(
    218218            PS_INTERPOLATE_BILINEAR,
    219             model->residuals->Ro, NULL, NULL, 0, 0.0, 0.0, 1, 0, 0.0);
     219            model->residuals->Ro, NULL, model->mask, 0, 0.0, 0.0, 1, 0, 0.0);
    220220        Rx = psImageInterpolateOptionsAlloc(
    221221            PS_INTERPOLATE_BILINEAR,
    … …  
    257257                float oy = yBin*(imageRow + 0.5 - yCenter) + yResidCenter;
    258258
     259                psU8 mflux = 0;
    259260                if (mode & PM_MODEL_OP_RES0) {
    260                     psU8 mflux = 0;
    261261                    double Fo = 0.0;
    262262                    psImageInterpolate (&Fo, NULL, &mflux, ox, oy, Ro);
    … …  
    265265                    }
    266266                }
    267                 if (mode & PM_MODEL_OP_RES1) {
    268                     psU8 mflux = 0;
     267                // skip Rx,Ry if Ro is masked
     268                if (!mflux && (mode & PM_MODEL_OP_RES1)) {
    269269                    double Fx = 0.0;
    270270                    double Fy = 0.0;

  • branches/eam_branch_20071023/psModules/src/objects/pmPSF_IO.c

    r15254 r15359  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2007-10-09 19:26:25 $
     8 *  @version $Revision: 1.27.2.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2007-10-23 20:54:53 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    9393}
    9494
    95 bool pmPSFmodelWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
    96 {
     95bool pmPSFmodelWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) {
    9796
    9897    pmFPA *fpa = file->fpa;

  • branches/eam_branch_20071023/psModules/src/objects/pmSource.h

    r15039 r15359  
    33 * @author EAM, IfA; GLG, MHPCC
    44 *
    5  * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
    6  * @date $Date: 2007-09-27 03:35:29 $
     5 * @version $Revision: 1.18.2.1 $ $Name: not supported by cvs2svn $
     6 * @date $Date: 2007-10-23 20:54:53 $
    77 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    88 */
    … …  
    8484    psRegion region;                    ///< area on image covered by selected pixels
    8585    float sky, skyErr;                  ///< The sky and its error at the center of the object
     86    pmSourceExtendedParameters *extpars; // extended source parameters
    8687};
     88
     89typedef struct {
     90  psVector *radius;
     91  psVector *flux;
     92} pmSourceRadialProfile;
     93
     94typedef struct {
     95  psVector *flux;
     96  psVector *fluxVar; // measured variance
     97  psVector *fluxErr; // formal error
     98} pmSourceAnnuli;
     99
     100typedef struct {
     101  float mag;
     102  float magErr;
     103  float rad;
     104  float radErr;
     105} pmSourceIsophotalValues;
     106
     107typedef struct {
     108  float mag;
     109  float magErr;
     110  float rad;
     111  float radErr;
     112} pmSourcePetrosianValues;
     113
     114typedef struct {
     115  float mag;
     116  float magErr;
     117  float rad;
     118  float radErr;
     119} pmSourceKronValues;
     120
     121typedef struct {
     122  pmSourceRadialProfile   *profile;
     123  pmSourceAnnuli          *annuli;
     124  pmSourceIsophotalValues *isophot;
     125  pmSourcePetrosianValues *petrosian;
     126  pmSourceKronValues      *kron;
     127} pmSourceExtendedParameters;
    87128
    88129/** pmPSFClump data structure

  • branches/eam_branch_20071023/psModules/src/objects/pmSourceIO.c

    r15227 r15359  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2007-10-05 22:46:34 $
     5 *  @version $Revision: 1.52.2.1 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2007-10-23 20:54:53 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    224224
    225225    bool status;
    226     char *exttype;
    227     char *dataname;
    228     char *headname;
    229226    pmHDU *hdu;
    230227    psMetadata *updates;
    231228    psMetadata *outhead;
     229
     230    char *exttype  = NULL;
     231    char *dataname = NULL;
     232    char *xsrcname = NULL;
     233    char *headname = NULL;
    232234
    233235    // XXX if sources is NULL, skip the cell or write out empty tables?
    … …  
    312314            }
    313315            dataname = pmFPAfileNameFromRule (rule, file, view);
     316
     317            if (XSRC_OUTPUT) {
     318              // EXTNAME for extended source data table
     319              rule = psMetadataLookupStr(&status, menu, "CMF.XSRC");
     320              if (!rule) {
     321                psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config");
     322                return false;
     323              }
     324              xsrcname = pmFPAfileNameFromRule (rule, file, view);
     325            }
    314326        }
    315327
    … …  
    375387            }
    376388            if (!strcmp (exttype, "PS1_DEV_1")) {
    377                 status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname);
    378             }
    379 
     389                status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname, xsrcname);
     390            }
    380391            if (!status) {
    381392                psError(PS_ERR_IO, false, "writing CMF data to %s with format %s\n", file->filename, exttype);

  • branches/eam_branch_20071023/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

    r15250 r15359  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2007-10-09 03:10:05 $
     5 *  @version $Revision: 1.4.2.1 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2007-10-23 20:54:53 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4747// XXX how do I generate the source tables which I need to send to PSPS?
    4848
    49 bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
     49bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, char *xsrcname)
    5050{
    5151
    … …  
    5757    psF32 xPos, yPos;
    5858    psF32 xErr, yErr;
     59
     60    // if we request XSRC output, add the XSRC name to this header
     61    if (xsrcname) {
     62      psMetadataAddStr (tableHeader, PS_LIST_TAIL, "EXTXSRC", PS_META_REPLACE, "name of XSRC table extension", xsrcname);
     63    }
    5964
    6065    // let's write these out in S/N order
    … …  
    144149        return false;
    145150    }
    146 
    147151    psFree (table);
     152
     153    if (xsrcname) {
     154      pmSourcesWriteXSRC_PS1_DEV_1 (file->fits, sources, xsrcname);
     155    }
     156
    148157    return true;
    149158}
    … …  
    241250    return sources;
    242251}
     252
     253bool pmSourcesWriteXSRC_PS1_DEV_1 (psFits *fits, psArray *sources, char *extname)
     254{
     255
     256    psArray *table;
     257    psMetadata *row;
     258    int i;
     259    psF32 *PAR, *dPAR;
     260    psEllipseAxes axes;
     261    psF32 xPos, yPos;
     262    psF32 xErr, yErr;
     263
     264    // create a header to hold the output data
     265    outhead = psMetadataAlloc ();
     266
     267    // write the links to the image header
     268    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);
     269
     270    // let's write these out in S/N order
     271    sources = psArraySort (sources, pmSourceSortBySN);
     272
     273    table = psArrayAllocEmpty (sources->n);
     274
     275    // we write out all sources, regardless of quality.  the source flags tell us the state
     276    for (i = 0; i < sources->n; i++) {
     277        pmSource *source = (pmSource *) sources->data[i];
     278        source->seq = i;
     279
     280        // skip source if it is not a ext sourc
     281
     282        // no difference between PSF and non-PSF model
     283        // XXX the PSF output should report the value for the psf, not the ext, model
     284        pmModel *model = source->modelEXT;
     285        if (model == NULL) continue;
     286
     287        PAR = model->params->data.F32;
     288        dPAR = model->dparams->data.F32;
     289        xPos = PAR[PM_PAR_XPOS];
     290        yPos = PAR[PM_PAR_YPOS];
     291        xErr = dPAR[PM_PAR_XPOS];
     292        yErr = dPAR[PM_PAR_YPOS];
     293
     294        axes = pmPSF_ModelToAxes (PAR, 20.0);
     295
     296        row = psMetadataAlloc ();
     297        // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
     298        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
     299        psMetadataAdd (row, PS_LIST_TAIL, "X_EXT",            PS_DATA_F32, "PSF x coordinate",                           xPos);
     300        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT",            PS_DATA_F32, "PSF y coordinate",                           yPos);
     301        psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr);
     302        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr);
     303        psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             PS_MIN (99.0, source->extMag));
     304        // XXX need to calculate psfMag, psfMagErr, extMag, extMagErr...
     305        psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        PS_MIN (99.0, source->errMag));
     306
     307        // XXX these should be major and minor, not 'x' and 'y'
     308        psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_X",      PS_DATA_F32, "PSF width in x coordinate",                  axes.major);
     309        psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_Y",      PS_DATA_F32, "PSF width in y coordinate",                  axes.minor);
     310        psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
     311
     312        // other values that I need to report:
     313        // R, Mag Petrosian, .....
     314
     315        psArrayAdd (table, 100, row);
     316        psFree (row);
     317    }
     318
     319    if (table->n == 0) {
     320        psFitsWriteBlank (fits, outhead, extname);
     321        psFree (table);
     322        return true;
     323    }
     324
     325    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     326    if (!psFitsWriteTable (fits, outhead, table, extname)) {
     327        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
     328        psFree(table);
     329        return false;
     330    }
     331    psFree (table);
     332
     333    return true;
     334}

  • branches/eam_branch_20071023/psastro/src/psastroAnalysis.c

    r15245 r15359  
    1616        return true;
    1717    }
    18 
    19     // XXX test retrun
    20     // return true;
    2118
    2219    // load the reference stars overlapping the data stars

  • branches/eam_branch_20071023/psastro/src/psastroDefineFiles.c

    r14212 r15359  
    44
    55    // these calls bind the I/O handle to the specified fpa
    6     pmFPAfile *output = NULL;
    76
    8     output = pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.OUTPUT");
     7    pmFPAfile *output = pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.OUTPUT");
    98    if (!output) {
    109        psError(PSASTRO_ERR_CONFIG, false, "Failed to build FPA from PSASTRO.INPUT");
    … …  
    1211    }
    1312    output->save = true;
     13
     14    // XXX not sure of this: should this be bind from args?
     15    pmFPAfile *refAstrom = pmFPAfileDefineFromConf (&status, config, "PSASTRO.REF.ASTROM");
     16    if (!refAstrom) {
     17        psError(PSASTRO_ERR_CONFIG, false, "Failed to load ref astrometry");
     18        return false;
     19    }
    1420
    1521# if (0)

  • branches/eam_branch_20071023/psastro/src/psastroMosaicAstrom.c

    r15260 r15359  
    2424    // compare chips with supplied mosaic model
    2525    // adjust significant outliers to match model
     26
     27    if (!psastroEnforceChips (fpa, recipe) {
     28        psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips");
     29        return false;
     30    }
    2631
    2732    // given the existing per-chip astrometry, determine matches between raw and ref stars

  • branches/eam_branch_20071023/psphot/src/psphotAnnuli.c

    r13983 r15359  
    33bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
    44
    5   psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
     5  assert (source->extpars);
     6  assert (source->extpars->profile);
     7  assert (source->extpars->profile->radius);
     8  assert (source->extpars->profile->flux);
     9
     10  psVector *radius = source->extpars->profile->radius;
     11  psVector *weight = source->extpars->profile->weight;
     12  psVector *flux = source->extpars->profile->flux;
     13
     14  // XXX how do I define the radii?  we can put a vector in the recipe...
     15  // radialBins defines the bounds or start and stop (we can skip some that way...
     16  psVector *radialBinsLower = psMetadataLookupVector (&status, recipe, RADIAL_ANNULAR_BINS_LOWER);
     17  psVector *radialBinsUpper = psMetadataLookupVector (&status, recipe, RADIAL_ANNULAR_BINS_UPPER);
     18  assert (radialBinsLower->n == radialBinsUpper->n);
     19
     20  psVector *fluxValues = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
     21  psVector *fluxSquare = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
     22  psVector *fluxWeight = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
     23  psVector *pixelCount = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
     24  psVectorInit (fluxValues, 0.0);
     25  psVectorInit (fluxSquare, 0.0);
     26  psVectorInit (fluxWeight, 0.0);
     27  psVectorInit (pixelCount, 0.0);
     28
     29  // XXX this code assumes the radii are in pixels.  convert from arcsec with plate scale
     30  // XXX assume the annulii above are not overlapping?  much faster...
     31  // XXX this might be must faster in the reverse order: loop over annulii and use disection to
     32  // skip to the start of the annulus.
     33  for (int i = 0; i < flux->n; i++) {
     34    for (int j = 0; j < radialBinsLower->n; j++) {
     35      if (radius->data.F32[i] < radialBinsLower->data.F32[j]) continue;
     36      if (radius->data.F32[i] > radialBinsUpper->data.F32[j]) continue;
     37      fluxValues->data.F32[j] += flux->data.F32[i];
     38      fluxSquare->data.F32[j] += PS_SQR(flux->data.F32[i]);
     39      fluxWeight->data.F32[j] += weight->data.F32[i];
     40      pixelCount->data.F32[j] += 1.0;
     41    }
     42  }
     43
     44  for (int j = 0; j < radialBinsLower->n; j++) {
     45    fluxValues->data.F32[j] /= pixelCount->data.F32[j];
     46    fluxSquare->data.F32[j] /= pixelCount->data.F32[j];
     47    fluxSquare->data.F32[j] -= PS_SQR(fluxValues->data.F32[j]);
     48  }
     49 
     50  source->extpars->annuli = pmSourceAnnuliAlloc ();
     51  source->extpars->annuli->flux    = fluxValues;
     52  source->extpars->annuli->fluxErr = fluxWeight;
     53  source->extpars->annuli->fluxVar = fluxSquare;
     54
    655  return true;
     56}
    757
    8 }

  • branches/eam_branch_20071023/psphot/src/psphotExtendedSources.c

    r13983 r15359  
    6666            return false;
    6767        } else {
     68          // XXX why am I caching the model?
    6869            pmSourceCacheModel (source, maskVal); // XXX put this in the source model function?
    6970            psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
    7071            Npsf ++;
     72        }
     73
     74        // all of these below require the radial profile
     75        // XXX push this as a test and call in each of the functions below?
     76        // XXX this constructs a pmSourceExtendedParameters element
     77        if (doPetrosian || doIsophotal || doAnnuli || doKron) {
     78          if (!psphotRadialProfile (source, recipe, maskVal)) {
     79            psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate radial profile");
     80            return false;
     81          }
    7182        }
    7283

  • branches/eam_branch_20071023/psphot/src/psphotIsophotal.c

    r13983 r15359  
    33bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
    44
    5   psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
     5  assert (source->extpars);
     6  assert (source->extpars->profile);
     7  assert (source->extpars->profile->radius);
     8  assert (source->extpars->profile->flux);
     9
     10  psVector *radius = source->extpars->profile->radius;
     11  psVector *flux = source->extpars->profile->flux;
     12
     13  // flux at which to measure isophotal parameters
     14  // XXX cache this?
     15  float ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX");
     16
     17  // find the first bin below the flux level and the last above the level
     18  // XXX can this be done faster with disection?
     19  // XXX do I need to worry about crazy outliers? 
     20  // XXX should i be smoothing or fitting the curve?
     21  firstBelow = -1;
     22  lastAbove = -1;
     23  for (int i = 0; i < flux->n; i++) {
     24    if (flux->data.F32[i] > ISOPHOT_FLUX) lastAbove = i;
     25    if ((firstBelow < 0) && (flux->data.F32[i] < ISOPHOT_FLUX)) firstBelow = i;
     26  }
     27
     28  // need to examine pixels in this vicinity
     29  float isophotalFluxFirst = 0;
     30  float isophotalFluxLast = 0;
     31  for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
     32    if (i <= firstBelow) {
     33      isophotFluxFirst += flux->data.F32[i];
     34    }
     35    if (i <= lastAbove) {
     36      isophotFluxLast += flux->data.F32[i];
     37    }
     38  }
     39  float isophotalFlux    = 0.5*(isophotFluxLast + isophotFluxFirst);
     40  float isophotalFluxErr = 0.5*fabs(isophotFluxLast - isophotFluxFirst);
     41
     42  float isophotalRad     = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]);
     43  float isophotalRadErr  = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]);
     44
     45  if (!source->extpars->isophot) {
     46    source->extpars->isophot = pmSourceIsophotalValuesAlloc ();
     47  }
     48 
     49  // these are uncalibrated: instrumental mags and pixel units
     50  source->extpars->isophot->mag    = -2.5*log10(isophotalFlux);
     51  source->extpars->isophot->magErr = isophotalFluxErr / isophotalFlux;
     52
     53  source->extpars->isophot->rad    = isophotalRad;
     54  source->extpars->isophot->radErr = isophotalRadErr;
     55
    656  return true;
    757
    858}
     59
     60# if (0)
     61  // XXX cache the tmpScalar
     62  psScalar fluxScalar;
     63  fluxScalar.type.type = PS_TYPE_F32;
     64  fluxScalar.data.F32  = ISOPHOT_FLUX;
     65
     66  // radius and flux must be sorted by radius
     67
     68  // XXX in general, flux decreases monotonically with radius.  however, since exceptions are
     69  // possible we need to do something to smooth or otherwise handle the flux variations
     70 
     71  // get the index of the flux-sorted vector
     72  psVector *fluxIndex = psVectorSortIndex (flux);
     73
     74  // XXX need to write psVectorIndexBinaryDisect, which operates on a
     75
     76  // use disection to get in the right vicinity
     77  binS = psVectorIndexBinaryDisect (&result, flux, fluxIndex, fluxScalar, PS_BISECT_FIRST);
     78  binE = psVectorIndexBinaryDisect (&result, flux, fluxIndex, fluxScalar, PS_BISECT_LAST);
     79
     80  // find the min and max radius bins in this range
     81
     82  // XXX do something to choose a more accurate radial bin
     83
     84# endif

  • branches/eam_branch_20071023/psphot/src/psphotMakeResiduals.c

    r14967 r15359  
    194194                resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
    195195                //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
     196
     197                // if value < NRESID_SIGMA*sigma, mask pixel in resid map
     198                if (resid->Ro->data.F32[oy][ox] < NRESID_SIGMA*fluxStats->sampleStdev) {
     199                  resid->mask->data.U8[oy][ox] = RESID_SIG_MASK;
     200                }
     201
    196202            } else {
    197203                assert (SPATIAL_ORDER == 1);
    … …  
    227233                resid->Rx->data.F32[oy][ox] = B->data.F64[1];
    228234                resid->Ry->data.F32[oy][ox] = B->data.F64[2];
     235
     236                dRo = sqrt(A->data.F32[0][0]);
     237                if (resid->Ro->data.F32[oy][ox] < NRESID_SIGMA*dRo) {
     238                  resid->mask->data.U8[oy][ox] = RESID_SIG_MASK;
     239                }
    229240                //resid->weight->data.F32[oy][ox] = XXX;
    230241            }

  • branches/eam_branch_20071023/psphot/src/psphotPetrosian.c

    r13983 r15359  
    33bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
    44
    5   psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
     5  assert (source->extpars);
     6  assert (source->extpars->profile);
     7  assert (source->extpars->profile->radius);
     8  assert (source->extpars->profile->flux);
     9
     10  psVector *radius = source->extpars->profile->radius;
     11  psVector *flux = source->extpars->profile->flux;
     12
     13  // flux at which to measure isophotal parameters
     14  // XXX cache this?
     15  float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
     16  float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
     17
     18  // first find flux at R0
     19  firstBelow = -1;
     20  lastAbove = -1;
     21  for (int i = 0; i < radius->n; i++) {
     22    if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i;
     23    if ((firstBelow < 0) && (flux->data.F32[i] < PETROSIAN_R0)) firstBelow = i;
     24  }
     25
     26  // average flux in this range
     27  fluxR0 = 0.0;
     28  fluxRn = 0;
     29  for (int i = PS_MIN(firstBelow, lastAbove); i <= PS_MAX(firstBelow, lastAbove); i++) {
     30    fluxR0 += flux->data.F32[i];
     31    fluxRn ++;
     32  }
     33  fluxR0 = fluxR0 / (float)(fluxRn);
     34 
     35  // target flux for petrosian radius
     36  fluxRP = fluxR0 * PETROSIAN_RF;
     37
     38  // find the first bin below the flux level and the last above the level
     39  // XXX can this be done faster with disection?
     40  // XXX do I need to worry about crazy outliers? 
     41  // XXX should i be smoothing or fitting the curve?
     42  firstBelow = -1;
     43  lastAbove = -1;
     44  for (int i = 0; i < flux->n; i++) {
     45    if (flux->data.F32[i] > fluxRP) lastAbove = i;
     46    if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i;
     47  }
     48
     49  // need to examine pixels in this vicinity
     50  float fluxFirst = 0;
     51  float fluxLast = 0;
     52  for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
     53    if (i <= firstBelow) {
     54      fluxFirst += flux->data.F32[i];
     55    }
     56    if (i <= lastAbove) {
     57      fluxLast += flux->data.F32[i];
     58    }
     59  }
     60  float flux    = 0.5*(fluxLast + fluxFirst);
     61  float fluxErr = 0.5*fabs(fluxLast - fluxFirst);
     62  // XXX need to use the weight appropriately here...
     63
     64  float rad     = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]);
     65  float radErr  = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]);
     66
     67  if (!source->extpars->petrosian) {
     68    source->extpars->petrosian = pmSourcePetrosianValuesAlloc ();
     69  }
     70 
     71  // these are uncalibrated: instrumental mags and pixel units
     72  source->extpars->petrosian->mag    = -2.5*log10(flux);
     73  source->extpars->petrosian->magErr = fluxErr / flux;
     74
     75  source->extpars->petrosian->rad    = rad;
     76  source->extpars->petrosian->radErr = radErr;
     77
    678  return true;
    779
Note: See TracChangeset for help on using the changeset viewer.