IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28013


Ignore:
Timestamp:
May 18, 2010, 4:11:53 PM (16 years ago)
Author:
eugene
Message:

merge changes from branches/eam_branches/psphot,psModules.20100506 (finish basic psphotStack)

Location:
trunk
Files:
63 edited
8 copied

Legend:

Unmodified
Added
Removed

  • trunk/psModules

  • trunk/psModules/src/objects/Makefile.am

    r27672 r28013  
    4343        pmSourceIO_CMF_PS1_V1.c \
    4444        pmSourceIO_CMF_PS1_V2.c \
     45        pmSourceIO_CMF_PS1_SV1.c \
    4546        pmSourceIO_CMF_PS1_DV1.c \
    4647        pmSourceIO_MatchedRefs.c \

  • trunk/psModules/src/objects/pmPhotObj.c

    r27657 r28013  
    5858    if (!object->sources) {
    5959        object->sources = psArrayAllocEmpty(1);
    60         object->x = source->peak->xf;
    61         object->y = source->peak->yf;
     60        object->x  = source->peak->xf;
     61        object->y  = source->peak->yf;
     62        object->SN = source->peak->SN;
     63    } else {
     64        object->SN = PS_MAX(object->SN, source->peak->SN);
    6265    }
    6366    psArrayAdd (object->sources, 1, source);
    6467    return true;
    6568}
     69
     70// sort by SN (descending)
     71int pmPhotObjSortBySN (const void **a, const void **b)
     72{
     73    pmPhotObj *objA = *(pmPhotObj **)a;
     74    pmPhotObj *objB = *(pmPhotObj **)b;
     75
     76    psF32 fA = objA->SN;
     77    psF32 fB = objB->SN;
     78
     79    psF32 diff = fA - fB;
     80    if (diff > FLT_EPSILON) return (-1);
     81    if (diff < FLT_EPSILON) return (+1);
     82    return (0);
     83}

  • trunk/psModules/src/objects/pmPhotObj.h

    r27657 r28013  
    3939    float x;
    4040    float y;
     41    float SN;                           // max of peak->SN for all matched sources
    4142} pmPhotObj;
    4243
     44pmPhotObj *pmPhotObjAlloc(void);
     45
    4346bool pmPhotObjAddSource(pmPhotObj *object, pmSource *source);
    44 pmPhotObj *pmPhotObjAlloc(void);
     47
     48int pmPhotObjSortBySN (const void **a, const void **b);
    4549
    4650/// @}

  • trunk/psModules/src/objects/pmSource.c

    r27657 r28013  
    5454    psFree(tmp->moments);
    5555    psFree(tmp->diffStats);
     56    psFree(tmp->radial);
    5657    psFree(tmp->blends);
    5758    psTrace("psModules.objects", 10, "---- end ----\n");
    … …  
    116117    source->extpars = NULL;
    117118    source->diffStats = NULL;
     119    source->radial = NULL;
    118120
    119121    source->region = psRegionSet(NAN, NAN, NAN, NAN);

  • trunk/psModules/src/objects/pmSource.h

    r27657 r28013  
    9191    pmSourceExtendedPars *extpars;      ///< extended source parameters
    9292    pmSourceDiffStats *diffStats;       ///< extra parameters for difference detections
     93    pmSourceRadialApertures *radial;    ///< radial flux in circular apertures
    9394    int imageID;
    9495};

  • trunk/psModules/src/objects/pmSourceExtendedPars.c

    r27818 r28013  
    6565}
    6666
     67// pmSourceRadialApertures carries the raw radial flux information, including angular bins
     68static void pmSourceRadialAperturesFree(pmSourceRadialApertures *radial)
     69{
     70    if (!radial) return;
     71    psFree(radial->flux);
     72    psFree(radial->fluxErr);
     73    psFree(radial->fill);
     74}
     75
     76pmSourceRadialApertures *pmSourceRadialAperturesAlloc()
     77{
     78    pmSourceRadialApertures *radial = (pmSourceRadialApertures *)psAlloc(sizeof(pmSourceRadialApertures));
     79    psMemSetDeallocator(radial, (psFreeFunc) pmSourceRadialAperturesFree);
     80
     81    radial->flux = NULL;
     82    radial->fluxErr = NULL;
     83    radial->fill = NULL;
     84    return radial;
     85}
     86
     87bool psMemCheckSourceRadialApertures(psPtr ptr)
     88{
     89    PS_ASSERT_PTR(ptr, false);
     90    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialAperturesFree);
     91}
     92
    6793// pmSourceEllipticalFlux carries the elliptical renormalized radial flux info
    6894static void pmSourceEllipticalFluxFree(pmSourceEllipticalFlux *flux)

  • trunk/psModules/src/objects/pmSourceExtendedPars.h

    r27818 r28013  
    2020    psVector *isophotalRadii;           // isophotal radius for the above angles
    2121} pmSourceRadialFlux;
     22
     23typedef struct {
     24    psVector *flux;                     // fluxes measured at above radii
     25    psVector *fluxErr;                  // formal error on the fluxes (sqrt\sum(variance))
     26    psVector *fill;                     // angles corresponding to above radial profiles
     27} pmSourceRadialApertures;
    2228
    2329typedef struct {
    … …  
    6268bool psMemCheckSourceRadialFlux(psPtr ptr);
    6369
     70pmSourceRadialApertures *pmSourceRadialAperturesAlloc();
     71bool psMemCheckSourceRadialApertures(psPtr ptr);
     72
    6473pmSourceEllipticalFlux *pmSourceEllipticalFluxAlloc();
    6574bool psMemCheckSourceEllipticalFlux(psPtr ptr);

  • trunk/psModules/src/objects/pmSourceIO.c

    r27818 r28013  
    539539                status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
    540540            }
     541            if (!strcmp (exttype, "PS1_SV1")) {
     542                status &= pmSourcesWrite_CMF_PS1_SV1 (file->fits, readout, sources, file->header, outhead, dataname, recipe);
     543            }
    541544            if (!strcmp (exttype, "PS1_DV1")) {
    542545                status &= pmSourcesWrite_CMF_PS1_DV1 (file->fits, readout, sources, file->header, outhead, dataname);
    … …  
    556559                    status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    557560                }
     561                if (!strcmp (exttype, "PS1_SV1")) {
     562                    status &= pmSourcesWrite_CMF_PS1_SV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
     563                }
    558564                if (!strcmp (exttype, "PS1_DV1")) {
    559565                    status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    … …  
    573579                    status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, readout, sources, xfitname);
    574580                }
     581                if (!strcmp (exttype, "PS1_SV1")) {
     582                    status &= pmSourcesWrite_CMF_PS1_SV1_XFIT (file->fits, readout, sources, xfitname);
     583                }
    575584                if (!strcmp (exttype, "PS1_DV1")) {
    576585                    status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, readout, sources, xfitname);
    … …  
    981990        // we need to find the corresponding table EXTNAME.
    982991        // first check the header
    983         char *extdata = psMetadataLookupStr (NULL, hdu->header, "EXTDATA");
     992        char *extdata = psMetadataLookupStr (&status, hdu->header, "EXTDATA");
    984993        if (extdata) {
    985994            // if EXTDATA is defined in the header, use that value for 'dataname'
    … …  
    10231032            if (!strcmp (exttype, "PS1_V2")) {
    10241033                sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header);
     1034            }
     1035            if (!strcmp (exttype, "PS1_SV1")) {
     1036                sources = pmSourcesRead_CMF_PS1_SV1 (file->fits, hdu->header);
    10251037            }
    10261038            if (!strcmp (exttype, "PS1_DV1")) {

  • trunk/psModules/src/objects/pmSourceIO.h

    r27818 r28013  
    4343bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname);
    4444
     45bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     46bool pmSourcesWrite_CMF_PS1_SV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     47bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname);
     48
    4549bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    4650bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    … …  
    5761psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header);
    5862psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header);
     63psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header);
    5964psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header);
    6065

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c

    r27818 r28013  
    175175        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    176176        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
    177         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",     PS_DATA_F32, "PSF fit instrumental magnitude",            source->psfFlux);
    178         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",       source->psfFluxErr);
     177        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfFlux);
     178        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfFluxErr);
    179179        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    180180        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);

  • trunk/psphot

  • trunk/psphot/doc/stack.txt

    r27848 r28013  
     1
     220100506:
     3
     4  we may load RAW and/or CNV images.
     5
     6  * options->convolve:
     7    * if only one of RAW or CNV exists, convolve it to match target PSF
     8    * if both exist, convolve desired one
     9   
     10    * if matching PSF exists, use it (unless told to re-measure)
     11   
     12    * output goes to OUT (which is then used by psphot analysis routines)
     13
     14  * detect
     15
     16    * if (RAW) ? RAW : OUT
     17   
     18  *   
     19
    120
    22120100503:

  • trunk/psphot/src/Makefile.am

    r27819 r28013  
    9595        psphotFitSourcesLinearStack.c \
    9696        psphotSourceMatch.c           \
     97        psphotStackMatchPSFs.c        \
     98        psphotStackMatchPSFsUtils.c   \
     99        psphotStackMatchPSFsPrepare.c \
     100        psphotStackOptions.c          \
     101        psphotStackObjects.c          \
     102        psphotStackPSF.c              \
    97103        psphotCleanup.c
    98104
    … …  
    159165        psphotModelWithPSF.c           \
    160166        psphotExtendedSourceAnalysis.c \
     167        psphotExtendedSourceAnalysisByObject.c \
    161168        psphotExtendedSourceFits.c     \
    162169        psphotKernelFromPSF.c          \
    … …  
    186193        psphotEllipticalProfile.c      \
    187194        psphotRadialBins.c             \
     195        psphotRadialApertures.c        \
     196        psphotRadialAperturesByObject.c \
    188197        psphotPetrosian.c              \
    189198        psphotPetrosianRadialBins.c    \

  • trunk/psphot/src/psphot.h

    r27819 r28013  
    3737bool            psphotReadoutMinimal(pmConfig *config, const pmFPAview *view);
    3838
    39 bool            psphotReadoutCleanup (pmConfig *config, const pmFPAview *view);
    40 bool            psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     39bool            psphotReadoutCleanup (pmConfig *config, const pmFPAview *view, const char *filerule);
     40bool            psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    4141
    4242bool            psphotDefineFiles (pmConfig *config, pmFPAfile *input);
    … …  
    5151
    5252// psphotReadout functions
    53 bool            psphotAddPhotcode (pmConfig *config, const pmFPAview *view);
     53bool            psphotAddPhotcode (pmConfig *config, const pmFPAview *view, const char *filerule);
    5454bool            psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    5555
    56 bool            psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view);
    57 bool            psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    58 
    59 bool            psphotModelBackground (pmConfig *config, const pmFPAview *view);
    60 bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filename, int index);
    61 
    62 bool            psphotSubtractBackground (pmConfig *config, const pmFPAview *view);
    63 bool            psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    64 
    65 bool            psphotFindDetections (pmConfig *config, const pmFPAview *view, bool firstPass);
    66 bool            psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool firstPass);
    67 
    68 bool            psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow);
    69 bool            psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool setWindow);
    70 
    71 bool            psphotDeblendSatstars (pmConfig *config, const pmFPAview *view);
    72 bool            psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
    73 
    74 bool            psphotBasicDeblend (pmConfig *config, const pmFPAview *view);
    75 bool            psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
    76 
    77 bool            psphotRoughClass (pmConfig *config, const pmFPAview *view);
    78 bool            psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     56bool            psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view, const char *filerule);
     57bool            psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     58
     59bool            psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
     60bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
     61
     62bool            psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
     63bool            psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     64
     65bool            psphotFindDetections (pmConfig *config, const pmFPAview *view, const char *filerule, bool firstPass);
     66bool            psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool firstPass);
     67
     68bool            psphotSourceStats (pmConfig *config, const pmFPAview *view, const char *filerule, bool setWindow);
     69bool            psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool setWindow);
     70
     71bool            psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule);
     72bool            psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
     73
     74bool            psphotBasicDeblend (pmConfig *config, const pmFPAview *view, const char *filerule);
     75bool            psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
     76
     77bool            psphotRoughClass (pmConfig *config, const pmFPAview *view, const char *filerule);
     78bool            psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    7979bool            psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *analysis, psMetadata *recipe, const bool havePSF);
    8080
    81 bool            psphotImageQuality (pmConfig *config, const pmFPAview *view);
    82 bool            psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    83 
    84 bool            psphotChoosePSF (pmConfig *config, const pmFPAview *view);
    85 bool            psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    86 
    87 bool            psphotGuessModels (pmConfig *config, const pmFPAview *view);
    88 bool            psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
    89 
    90 bool            psphotMergeSources (pmConfig *config, const pmFPAview *view);
    91 bool            psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
    92 
    93 bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final);
     81bool            psphotImageQuality (pmConfig *config, const pmFPAview *view, const char *filerule);
     82bool            psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     83
     84bool            psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule);
     85bool            psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     86
     87bool            psphotGuessModels (pmConfig *config, const pmFPAview *view, const char *filerule);
     88bool            psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
     89
     90bool            psphotMergeSources (pmConfig *config, const pmFPAview *view, const char *filerule);
     91bool            psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
     92
     93bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
    9494bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
    9595
    96 bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize);
    97 bool            psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool getPSFsize);
    98 
    99 bool            psphotBlendFit (pmConfig *config, const pmFPAview *view);
    100 bool            psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     96bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);
     97bool            psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool getPSFsize);
     98
     99bool            psphotBlendFit (pmConfig *config, const pmFPAview *view, const char *filerule);
     100bool            psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    101101bool            psphotBlendFit_Threaded (psThreadJob *job);
    102102
    103 bool            psphotReplaceAllSources (pmConfig *config, const pmFPAview *view);
    104 bool            psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    105 
    106 bool            psphotAddNoise (pmConfig *config, const pmFPAview *view);
    107 bool            psphotSubNoise (pmConfig *config, const pmFPAview *view);
    108 bool            psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add);
    109 bool            psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add);
    110 
    111 bool            psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view);
    112 bool            psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    113 
    114 bool            psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view);
    115 bool            psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    116 
    117 bool            psphotApResid (pmConfig *config, const pmFPAview *view);
    118 bool            psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    119 
    120 bool            psphotMagnitudes (pmConfig *config, const pmFPAview *view);
     103bool            psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule);
     104bool            psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     105
     106bool            psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule);
     107bool            psphotSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule);
     108bool            psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule, bool add);
     109bool            psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add);
     110
     111bool            psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule);
     112bool            psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     113
     114bool            psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule);
     115bool            psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     116
     117bool            psphotApResid (pmConfig *config, const pmFPAview *view, const char *filerule);
     118bool            psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     119
     120bool            psphotMagnitudes (pmConfig *config, const pmFPAview *view, const char *filerule);
    121121bool            psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
    122122bool            psphotMagnitudes_Threaded (psThreadJob *job);
    123123
    124 bool            psphotEfficiency (pmConfig *config, const pmFPAview *view);
    125 bool            psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     124bool            psphotEfficiency (pmConfig *config, const pmFPAview *view, const char *filerule);
     125bool            psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    126126
    127127bool            psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources);
    128128bool            psphotPSFWeights_Threaded (psThreadJob *job);
    129129
    130 bool            psphotSkyReplace (pmConfig *config, const pmFPAview *view);
    131 bool            psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
    132 
    133 bool            psphotSourceFreePixels (pmConfig *config, const pmFPAview *view);
    134 bool            psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index);
     130bool            psphotSkyReplace (pmConfig *config, const pmFPAview *view, const char *filerule);
     131bool            psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
     132
     133bool            psphotSourceFreePixels (pmConfig *config, const pmFPAview *view, const char *filerule);
     134bool            psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    135135
    136136// in psphotSourceStats.c:
    … …  
    147147
    148148// in psphotMergeSources.c:
    149 bool            psphotLoadExtSources (pmConfig *config, const pmFPAview *view);
     149bool            psphotLoadExtSources (pmConfig *config, const pmFPAview *view, const char *filerule);
    150150psArray        *psphotLoadPSFSources (pmConfig *config, const pmFPAview *view);
    151 bool            psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view);
    152 bool            psphotCheckExtSources (pmConfig *config, const pmFPAview *view);
     151bool            psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view, const char *filerule);
     152bool            psphotCheckExtSources (pmConfig *config, const pmFPAview *view, const char *filerule);
    153153
    154154// generate the detection structure for the supplied array of sources
    155 bool            psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, psArray *sources);
     155bool            psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *sources);
    156156
    157157// generate the detection structure for the supplied array of sources
    … …  
    349349bool psphotStackImageLoop (pmConfig *config);
    350350bool psphotStackReadout (pmConfig *config, const pmFPAview *view);
    351 bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view);
     351bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleCnv);
    352352bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration
    353353                                     const pmFPAview *view,
     354                                     const char *filename,
    354355                                     pmReadout **chiReadout,
    355                                      char *filename,
    356356                                     int index);
    357357
    358 bool psphotStackRemoveChisqFromInputs (pmConfig *config);
     358bool psphotStackRemoveChisqFromInputs (pmConfig *config, const char *filerule);
    359359bool pmFPAfileRemoveSingle(psMetadata *files, const char *name, int num);
    360360
    361 psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view);
    362 bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index);
     361psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view, const char *filerule);
     362bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    363363bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS);
    364364
    … …  
    367367int pmPhotObjSortByX (const void **a, const void **b);
    368368
     369typedef enum {
     370    PSPHOT_CNV_SRC_NONE,
     371    PSPHOT_CNV_SRC_AUTO,
     372    PSPHOT_CNV_SRC_CNV,
     373    PSPHOT_CNV_SRC_RAW,
     374} psphotStackConvolveSource;
     375
     376/// Options for stacking process
     377typedef struct {
     378    // Setup
     379   
     380    int numCols;                            // size of image (X)
     381    int numRows;                            // size of image (Y)
     382
     383    int num;                            // Number of inputs
     384    bool convolve;                      // Convolve images?
     385    psphotStackConvolveSource convolveSource;
     386    // psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images
     387
     388    // bool matchZPs;                      // Adjust relative fluxes based on transparency analysis?
     389    // bool photometry;                    // Perform photometry?
     390    // psMetadata *stats;                  // Statistics for output
     391    // FILE *statsFile;                    // File to which to write statistics
     392    // psArray *origImages, *origMasks, *origVariances; // Filenames of the original images
     393    // psArray *origCovars;                // Original covariances matrices
     394    // int quality;                        // Bad data quality flag
     395
     396    // Prepare
     397    pmPSF *psf;                         // Target PSF
     398    psVector *inputSeeing;              // Input seeing FWHMs
     399    psVector *inputMask;                // Mask for inputs
     400
     401    float targetSeeing;                 // Target seeing FWHM
     402    psArray *sourceLists;               // Individual lists of sources for matching
     403    psVector *norm;                     // Normalisation for each image
     404    psArray *psfs;
     405
     406    // psVector *exposures;                // Exposure times
     407    // float sumExposure;                  // Sum of exposure times
     408    // float zp;                           // Zero point for output
     409    // psVector *inputMask;                // Mask for inputs
     410    // psArray *sources;                   // Matched sources
     411
     412    // Convolve
     413    psArray *kernels;                   // PSF-matching kernels --- required in the stacking
     414    psArray *regions;                   // PSF-matching regions --- required in the stacking
     415    psVector *matchChi2;                // chi^2 for stamps from matching
     416    psVector *weightings;               // Combination weightings for images (1/noise^2)
     417    // psArray *cells;                     // Cells for convolved images --- a handle for reading again
     418    // int numCols, numRows;               // Size of image
     419    // psArray *convCovars;                // Convolved covariance matrices
     420
     421    // Combine initial
     422    // pmReadout *outRO;                   // Output readout
     423    // pmReadout *expRO;                   // Exposure readout
     424    // psArray *inspect;                   // Array of arrays of pixels to inspect
     425
     426    // Rejection
     427    // psArray *rejected;                  // Rejected pixels
     428} psphotStackOptions;
     429
     430/*** psphotStackMatchPSF prototypes ***/
     431bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view);
     432bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
     433bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
     434
     435// psphotStackMatchPSFsUtils
     436psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
     437pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
     438bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index);
     439bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index);
     440bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index);
     441bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index);
     442bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
     443bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
     444bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index);
     445
     446pmPSF *psphotStackPSF(const pmConfig *config, int numCols, int numRows, const psArray *psfs, const psVector *inputMask);
     447
     448psphotStackOptions *psphotStackOptionsAlloc (int num);
     449psphotStackConvolveSource psphotStackConvolveSourceFromString (const char *string);
     450pmFPAfile *psphotStackGetConvolveSource (pmConfig *config, psphotStackOptions *options, int index);
     451
     452bool psphotCopySources (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
     453bool psphotCopySourcesReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
     454
     455bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule);
     456bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     457bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax);
     458
     459bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
     460bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
     461
     462bool psphotStackObjectsUnifyPosition (psArray *objects);
     463
    369464#endif

  • trunk/psphot/src/psphotAddNoise.c

    r26894 r28013  
    11# include "psphotInternal.h"
    22
    3 bool psphotAddNoise (pmConfig *config, const pmFPAview *view) {
    4     return psphotAddOrSubNoise (config, view, true);
     3bool psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule) {
     4    return psphotAddOrSubNoise (config, view, filerule, true);
    55}
    66
    7 bool psphotSubNoise (pmConfig *config, const pmFPAview *view) {
    8     return psphotAddOrSubNoise (config, view, false);
     7bool psphotSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule) {
     8    return psphotAddOrSubNoise (config, view, filerule, false);
    99}
    1010
    1111// for now, let's store the detections on the readout->analysis for each readout
    12 bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add)
     12bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule, bool add)
    1313{
    1414    bool status = true;
    … …  
    2323    // loop over the available readouts
    2424    for (int i = 0; i < num; i++) {
    25         if (!psphotAddOrSubNoiseReadout (config, view, "PSPHOT.INPUT", i, recipe, add)) {
    26             psError (PSPHOT_ERR_CONFIG, false, "failed on to modify noise for PSPHOT.INPUT entry %d", i);
     25        if (!psphotAddOrSubNoiseReadout (config, view, filerule, i, recipe, add)) {
     26            psError (PSPHOT_ERR_CONFIG, false, "failed on to modify noise for %s entry %d", filerule, i);
    2727            return false;
    2828        }
    … …  
    3131}
    3232
    33 bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add) {
     33bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add) {
    3434
    3535    bool status = false;
    … …  
    3939
    4040    // find the currently selected readout
    41     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     41    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    4242    psAssert (file, "missing file?");
    4343

  • trunk/psphot/src/psphotApResid.c

    r27657 r28013  
    55
    66// for now, let's store the detections on the readout->analysis for each readout
    7 bool psphotApResid (pmConfig *config, const pmFPAview *view)
     7bool psphotApResid (pmConfig *config, const pmFPAview *view, const char *filerule)
    88{
    99    bool status = true;
    … …  
    2323    for (int i = 0; i < num; i++) {
    2424        if (i == chisqNum) continue; // skip chisq image
    25         if (!psphotApResidReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    26             psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for PSPHOT.INPUT entry %d", i);
     25        if (!psphotApResidReadout (config, view, filerule, i, recipe)) {
     26            psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for %s entry %d", filerule, i);
    2727            return false;
    2828        }
    … …  
    3131}
    3232
    33 bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
     33bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe)
    3434{
    3535    int Nfail = 0;
    … …  
    4343
    4444    // find the currently selected readout
    45     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     45    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    4646    psAssert (file, "missing file?");
    4747

  • trunk/psphot/src/psphotBasicDeblend.c

    r26894 r28013  
    22
    33// for now, let's store the detections on the readout->analysis for each readout
    4 bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view)
     4bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view, const char *filerule)
    55{
    66    bool status = true;
    … …  
    1111    // loop over the available readouts
    1212    for (int i = 0; i < num; i++) {
    13         if (!psphotBasicDeblendReadout (config, view, "PSPHOT.INPUT", i)) {
    14             psError (PSPHOT_ERR_CONFIG, false, "failed on basic deblend analysis for PSPHOT.INPUT entry %d", i);
     13        if (!psphotBasicDeblendReadout (config, view, filerule, i)) {
     14            psError (PSPHOT_ERR_CONFIG, false, "failed on basic deblend analysis for %s entry %d", filerule, i);
    1515            return false;
    1616        }
    … …  
    1919}
    2020
    21 bool psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filename, int fileIndex) {
     21bool psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) {
    2222
    2323    int N;
    … …  
    3131
    3232    // find the currently selected readout
    33     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, fileIndex); // File of interest
     33    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // File of interest
    3434    psAssert (file, "missing file?");
    3535

  • trunk/psphot/src/psphotBlendFit.c

    r26894 r28013  
    22
    33// for now, let's store the detections on the readout->analysis for each readout
    4 bool psphotBlendFit (pmConfig *config, const pmFPAview *view)
     4bool psphotBlendFit (pmConfig *config, const pmFPAview *view, const char *filerule)
    55{
    66    bool status = true;
    … …  
    1515    // loop over the available readouts
    1616    for (int i = 0; i < num; i++) {
    17         if (!psphotBlendFitReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    18             psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (non-linear) for PSPHOT.INPUT entry %d", i);
     17        if (!psphotBlendFitReadout (config, view, filerule, i, recipe)) {
     18            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (non-linear) for %s entry %d", filerule, i);
    1919            return false;
    2020        }
    … …  
    2424
    2525// XXX I don't like this name
    26 bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     26bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
    2727
    2828    int Nfit = 0;
    … …  
    3535
    3636    // find the currently selected readout
    37     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     37    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    3838    psAssert (file, "missing file?");
    3939

  • trunk/psphot/src/psphotChoosePSF.c

    r27657 r28013  
    22
    33// generate a PSF model for inputs without PSF models already loaded
    4 bool psphotChoosePSF (pmConfig *config, const pmFPAview *view)
     4bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule)
    55{
    66    bool status = true;
    … …  
    2020    for (int i = 0; i < num; i++) {
    2121        if (i == chisqNum) continue; // skip chisq image
    22         if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    23             psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for PSPHOT.INPUT entry %d", i);
     22        if (!psphotChoosePSFReadout (config, view, filerule, i, recipe)) {
     23            psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for %s entry %d", filerule, i);
    2424            return false;
    2525        }
    … …  
    2929
    3030// try PSF models and select best option
    31 bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     31bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
    3232
    3333    bool status;
    … …  
    3636
    3737    // find the currently selected readout
    38     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     38    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    3939    psAssert (file, "missing file?");
    4040

  • trunk/psphot/src/psphotCleanup.c

    r27833 r28013  
    1919    pmConceptsDone ();
    2020    pmConfigDone ();
     21    psLibFinalize();
    2122    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot");
    2223    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot");

  • trunk/psphot/src/psphotDeblendSatstars.c

    r26894 r28013  
    22
    33// for now, let's store the detections on the readout->analysis for each readout
    4 bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view)
     4bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule)
    55{
    66    bool status = true;
    … …  
    1111    // loop over the available readouts
    1212    for (int i = 0; i < num; i++) {
    13         if (!psphotDeblendSatstarsReadout (config, view, "PSPHOT.INPUT", i)) {
    14             psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     13        if (!psphotDeblendSatstarsReadout (config, view, filerule, i)) {
     14            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for %s entry %d", filerule, i);
    1515            return false;
    1616        }
    … …  
    1919}
    2020
    21 bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int fileIndex) {
     21bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) {
    2222
    2323    int N;
    … …  
    3131
    3232    // find the currently selected readout
    33     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, fileIndex); // File of interest
     33    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // File of interest
    3434    psAssert (file, "missing file?");
    3535

  • trunk/psphot/src/psphotEfficiency.c

    r28006 r28013  
    156156}
    157157
    158 bool psphotEfficiency (pmConfig *config, const pmFPAview *view)
     158bool psphotEfficiency (pmConfig *config, const pmFPAview *view, const char *filerule)
    159159{
    160160    bool status = true;
    … …  
    173173    // loop over the available readouts
    174174    for (int i = 0; i < num; i++) {
    175         if (i == chisqNum) continue; // skip chisq image
    176         if (!psphotEfficiencyReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    177             psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for PSPHOT.INPUT entry %d", i);
     175        if (i == chisqNum) continue; // skip chisq image
     176        if (!psphotEfficiencyReadout (config, view, filerule, i, recipe)) {
     177            psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for %s entry %d", filerule, i);
    178178            return false;
    179179        }
    … …  
    183183
    184184// Determine detection efficiency
    185 bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
     185bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe)
    186186{
    187187    bool status = true;
    … …  
    190190
    191191    // find the currently selected readout
    192     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     192    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    193193    psAssert (file, "missing file?");
    194194

  • trunk/psphot/src/psphotErrorCodes.dat

    r27002 r28013  
    1111APERTURE                Problem with aperture photometry
    1212SKY                     Problem in sky determination
     13IO                      Problem in data I/O
    1314# these errors correspond to standard exit conditions
    1415ARGUMENTS               Incorrect arguments

  • trunk/psphot/src/psphotExtendedSourceAnalysis.c

    r27819 r28013  
    11# include "psphotInternal.h"
    22
     3// ?? these cannot happen here --> we would need to do this in psphotExtendedSourceAnalysis
     4// XXX option to choose a consistent position
     5// XXX option to choose a consistent elliptical contour
     6// XXX SDSS uses the r-band petrosian radius to measure petrosian fluxes in all bands
     7// XXX consistent choice of extendedness...
     8
    39// for now, let's store the detections on the readout->analysis for each readout
    4 bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view)
     10bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule)
    511{
    612    bool status = true;
    … …  
    2127    // loop over the available readouts
    2228    for (int i = 0; i < num; i++) {
    23         if (!psphotExtendedSourceAnalysisReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    24             psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for PSPHOT.INPUT entry %d", i);
     29        if (!psphotExtendedSourceAnalysisReadout (config, view, filerule, i, recipe)) {
     30            psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
    2531            return false;
    2632        }
    … …  
    3036
    3137// aperture-like measurements for extended sources
    32 bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     38bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
    3339
    3440    bool status;
    … …  
    4248
    4349    // find the currently selected readout
    44     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     50    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    4551    psAssert (file, "missing file?");
    4652

  • trunk/psphot/src/psphotExtendedSourceFits.c

    r26894 r28013  
    22
    33// for now, let's store the detections on the readout->analysis for each readout
    4 bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view)
     4bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule)
    55{
    66    bool status = true;
    … …  
    2121    // loop over the available readouts
    2222    for (int i = 0; i < num; i++) {
    23         if (!psphotExtendedSourceFitsReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    24             psError (PSPHOT_ERR_CONFIG, false, "failed on to fit extended sources for PSPHOT.INPUT entry %d", i);
     23        if (!psphotExtendedSourceFitsReadout (config, view, filerule, i, recipe)) {
     24            psError (PSPHOT_ERR_CONFIG, false, "failed on to fit extended sources for %s entry %d", filerule, i);
    2525            return false;
    2626        }
    … …  
    3030
    3131// non-linear model fitting for extended sources
    32 bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     32bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
    3333
    3434    bool status;
    … …  
    4141
    4242    // find the currently selected readout
    43     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     43    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    4444    psAssert (file, "missing file?");
    4545

  • trunk/psphot/src/psphotFindDetections.c

    r26894 r28013  
    44// peaks and new footprints.  any old peaks are saved on oldPeaks.  the resulting footprint set
    55// contains all footprints (old and new)
    6 bool psphotFindDetections (pmConfig *config, const pmFPAview *view, bool firstPass)
     6bool psphotFindDetections (pmConfig *config, const pmFPAview *view, const char *filerule, bool firstPass)
    77{
    88    bool status = true;
    … …  
    1717    // loop over the available readouts
    1818    for (int i = 0; i < num; i++) {
    19         if (!psphotFindDetectionsReadout (config, view, "PSPHOT.INPUT", i, recipe, firstPass)) {
    20             psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
     19        if (!psphotFindDetectionsReadout (config, view, filerule, i, recipe, firstPass)) {
     20            psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for %s entry %d", filerule, i);
    2121            return false;
    2222        }
    … …  
    2626
    2727// smooth the image, search for peaks, optionally define footprints based on the peaks
    28 bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool firstPass) {
     28bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool firstPass) {
    2929
    3030    bool status;
    … …  
    3434
    3535    // find the currently selected readout
    36     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     36    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    3737    psAssert (file, "missing file?");
    3838

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r27532 r28013  
    1313
    1414// for now, let's store the detections on the readout->analysis for each readout
    15 bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final)
     15bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final)
    1616{
    1717    bool status = true;
    … …  
    2828
    2929        // find the currently selected readout
    30         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest
     30        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
    3131        psAssert (file, "missing file?");
    3232
    … …  
    4444
    4545        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) {
    46             psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i);
     46            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
    4747            return false;
    4848        }

  • trunk/psphot/src/psphotFitSourcesLinearStack.c

    r27657 r28013  
    2828
    2929    // analysis is done in spatial order (to speed up overlap search)
    30     // sort by first element in each source list
    3130    objects = psArraySort (objects, pmPhotObjSortByX);
    3231

  • trunk/psphot/src/psphotForcedReadout.c

    r27314 r28013  
    2020
    2121    // set the photcode for this image
    22     if (!psphotAddPhotcode (config, view)) {
     22    if (!psphotAddPhotcode (config, view, "PSPHOT.INPUT")) {
    2323        psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
    2424        return false;
    … …  
    3030
    3131    // Generate the mask and weight images, including the user-defined analysis region of interest
    32     psphotSetMaskAndVariance (config, view);
     32    psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT");
    3333    if (!strcasecmp (breakPt, "NOTHING")) {
    34         return psphotReadoutCleanup(config, view);
     34        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    3535    }
    3636
    3737    // generate a background model (median, smoothed image)
    38     if (!psphotModelBackground (config, view)) {
    39         return psphotReadoutCleanup (config, view);
     38    if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) {
     39        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4040    }
    41     if (!psphotSubtractBackground (config, view)) {
    42         return psphotReadoutCleanup (config, view);
     41    if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) {
     42        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4343    }
    4444    if (!strcasecmp (breakPt, "BACKMDL")) {
    45         return psphotReadoutCleanup (config, view);
     45        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4646    }
    4747
    … …  
    4949        // this only happens if we had a programming error in psphotLoadPSF
    5050        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
    51         return psphotReadoutCleanup (config, view);
     51        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    5252    }
    5353
    5454    // include externally-supplied sources
    55     psphotLoadExtSources (config, view);
     55    psphotLoadExtSources (config, view, "PSPHOT.INPUT");
    5656
    5757    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    58     psphotGuessModels (config, view);
     58    psphotGuessModels (config, view, "PSPHOT.INPUT");
    5959
    6060    // merge the newly selected sources into the existing list
    6161    // NOTE: merge OLD and NEW
    62     psphotMergeSources (config, view);
     62    psphotMergeSources (config, view, "PSPHOT.INPUT");
    6363
    6464    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    65     psphotFitSourcesLinear (config, view, false);
     65    psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", false);
    6666
    6767    // identify CRs and extended sources
    … …  
    7171
    7272    // calculate source magnitudes
    73     psphotMagnitudes(config, view);
     73    psphotMagnitudes(config, view, "PSPHOT.INPUT");
    7474
    7575    // XXX do I want to do this?
    … …  
    8080
    8181    // replace background in residual image
    82     psphotSkyReplace (config, view);
     82    psphotSkyReplace (config, view, "PSPHOT.INPUT");
    8383
    8484    // drop the references to the image pixels held by each source
    85     psphotSourceFreePixels (config, view);
     85    psphotSourceFreePixels (config, view, "PSPHOT.INPUT");
    8686
    8787    // create the exported-metadata and free local data
    88     return psphotReadoutCleanup(config, view);
     88    return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    8989}

  • trunk/psphot/src/psphotGuessModels.c

    r27657 r28013  
    88
    99// for now, let's store the detections on the readout->analysis for each readout
    10 bool psphotGuessModels (pmConfig *config, const pmFPAview *view)
     10bool psphotGuessModels (pmConfig *config, const pmFPAview *view, const char *filerule)
    1111{
    1212    bool status = true;
    … …  
    2222    for (int i = 0; i < num; i++) {
    2323        if (i == chisqNum) continue; // skip chisq image
    24         if (!psphotGuessModelsReadout (config, view, "PSPHOT.INPUT", i)) {
    25             psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for PSPHOT.INPUT entry %d", i);
     24        if (!psphotGuessModelsReadout (config, view, filerule, i)) {
     25            psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for %s entry %d", filerule, i);
    2626            return false;
    2727        }
    … …  
    3131
    3232// construct an initial PSF model for each object (new sources only)
    33 bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
     33bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
    3434
    3535    bool status;
    … …  
    3838
    3939    // find the currently selected readout
    40     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     40    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    4141    psAssert (file, "missing file?");
    4242

  • trunk/psphot/src/psphotImageQuality.c

    r27657 r28013  
    55
    66// for now, let's store the detections on the readout->analysis for each readout
    7 bool psphotImageQuality (pmConfig *config, const pmFPAview *view)
     7bool psphotImageQuality (pmConfig *config, const pmFPAview *view, const char *filerule)
    88{
    99    bool status = true;
    … …  
    2323    for (int i = 0; i < num; i++) {
    2424        if (i == chisqNum) continue; // skip chisq image
    25         if (!psphotImageQualityReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    26             psError (PSPHOT_ERR_CONFIG, false, "failed on to measure image quality for PSPHOT.INPUT entry %d", i);
     25        if (!psphotImageQualityReadout (config, view, filerule, i, recipe)) {
     26            psError (PSPHOT_ERR_CONFIG, false, "failed on to measure image quality for %s entry %d", filerule, i);
    2727            return false;
    2828        }
    … …  
    3232
    3333// selecting the 'good' stars (likely to be psf stars), measure the M_cn, M_sn terms for n = 2,3,4
    34 bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
     34bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe)
    3535{
    3636    bool status = true;
    3737
    3838    // find the currently selected readout
    39     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     39    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    4040    psAssert (file, "missing file?");
    4141

  • trunk/psphot/src/psphotMagnitudes.c

    r27657 r28013  
    11# include "psphotInternal.h"
    22
    3 bool psphotMagnitudes (pmConfig *config, const pmFPAview *view)
     3bool psphotMagnitudes (pmConfig *config, const pmFPAview *view, const char *filerule)
    44{
    55    bool status = true;
    … …  
    2121
    2222        // find the currently selected readout
    23         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest
     23        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
    2424        psAssert (file, "missing file?");
    2525
    … …  
    3737
    3838        if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) {
    39             psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for PSPHOT.INPUT entry %d", i);
     39            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
    4040            return false;
    4141        }

  • trunk/psphot/src/psphotMakePSFReadout.c

    r26894 r28013  
    1919
    2020    // set the photcode for this image
    21     if (!psphotAddPhotcode (config, view)) {
     21    if (!psphotAddPhotcode (config, view, "PSPHOT.INPUT")) {
    2222        psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
    2323        return false;
    … …  
    2929
    3030    // Generate the mask and weight images, including the user-defined analysis region of interest
    31     psphotSetMaskAndVariance (config, view);
     31    psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT");
    3232    if (!strcasecmp (breakPt, "NOTHING")) {
    33         return psphotReadoutCleanup(config, view);
     33        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    3434    }
    3535
    3636    // generate a background model (median, smoothed image)
    37     if (!psphotModelBackground (config, view)) {
    38         return psphotReadoutCleanup (config, view);
     37    if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) {
     38        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    3939    }
    40     if (!psphotSubtractBackground (config, view)) {
    41         return psphotReadoutCleanup (config, view);
     40    if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) {
     41        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4242    }
    4343    if (!strcasecmp (breakPt, "BACKMDL")) {
    44         return psphotReadoutCleanup (config, view);
     44        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4545    }
    4646
    47     psphotLoadExtSources (config, view);
     47    psphotLoadExtSources (config, view, "PSPHOT.INPUT");
    4848
    4949    // If sources have been supplied, then these should be used to measure the PSF include
    … …  
    5353    // a text file have no valid SN values, but psphotChoosePSF needs to select the top
    5454    // PSF_MAX_NSTARS to generate the PSF.
    55     if (!psphotCheckExtSources (config, view)) {
     55    if (!psphotCheckExtSources (config, view, "PSPHOT.INPUT")) {
    5656        psLogMsg ("psphot", 3, "failure to select possible PSF sources (external or internal)");
    57         return psphotReadoutCleanup (config, view);
     57        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    5858    }
    5959
    6060    // Use bright stellar objects to measure PSF. If we do not have enough stars to generate
    6161    // the PSF, build one from the SEEING guess and model class
    62     if (!psphotChoosePSF (config, view)) {
     62    if (!psphotChoosePSF (config, view, "PSPHOT.INPUT")) {
    6363        psLogMsg ("psphot", 3, "failure to construct a psf model");
    64         return psphotReadoutCleanup (config, view);
     64        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    6565    }
    6666
    6767    // measure aperture photometry corrections
    6868# if 0
    69     if (!psphotApResid (config, view)) {
     69    if (!psphotApResid (config, view, "PSPHOT.INPUT")) {
    7070        psLogMsg ("psphot", 3, "failed on psphotApResid");
    71         return psphotReadoutCleanup (config, view);
     71        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    7272    }
    7373# endif
    7474
    75     return psphotReadoutCleanup (config, view);
     75    return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    7676}

  • trunk/psphot/src/psphotMaskReadout.c

    r26894 r28013  
    11# include "psphotInternal.h"
    22
    3 bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view) {
     3bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view, const char *filerule) {
    44
    55    bool status = false;
    … …  
    1616
    1717        // Generate the mask and weight images, including the user-defined analysis region of interest
    18         if (!psphotSetMaskAndVarianceReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     18        if (!psphotSetMaskAndVarianceReadout (config, view, filerule, i, recipe)) {
    1919            psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for PSPHOT.INPUT entry %d", i);
    2020            return false;
    … …  
    2525
    2626// generate mask and variance if not defined, additional mask for restricted subregion
    27 bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     27bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
    2828
    2929    bool status;
    3030
    31     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
     31    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    3232    psAssert (file, "missing file?");
    3333

  • trunk/psphot/src/psphotMergeSources.c

    r27657 r28013  
    66
    77// for now, let's store the detections on the readout->analysis for each readout
    8 bool psphotMergeSources (pmConfig *config, const pmFPAview *view)
     8bool psphotMergeSources (pmConfig *config, const pmFPAview *view, const char *filerule)
    99{
    1010    bool status = true;
    … …  
    1515    // loop over the available readouts
    1616    for (int i = 0; i < num; i++) {
    17         if (!psphotMergeSourcesReadout (config, view, "PSPHOT.INPUT", i)) {
    18             psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i);
     17        if (!psphotMergeSourcesReadout (config, view, filerule, i)) {
     18            psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for %s entry %d", filerule, i);
    1919            return false;
    2020        }
    … …  
    2424
    2525// add newly selected sources to the existing list of sources
    26 bool psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
     26bool psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
    2727
    2828    bool status;
    2929
    3030    // find the currently selected readout
    31     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     31    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    3232    psAssert (file, "missing file?");
    3333
    … …  
    7171// only expect a single entry for PSPHOT.INPUT.CMF and PSPHOT.SOURCES.TEXT, so we can only
    7272// associate input sources with a single entry for PSPHOT.INPUT
    73 bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view) {
     73bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view, const char *filerule) {
    7474
    7575    bool status;
    … …  
    7979
    8080    // find the currently selected readout
    81     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
     81    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    8282    psAssert (file, "missing file?");
    8383
    … …  
    130130    // load data from input TXT file:
    131131    {
    132         pmChip *chipTXT = pmFPAfileThisChip (config->files, view, "PSPHOT.INPUT");
     132        pmChip *chipTXT = pmFPAfileThisChip (config->files, view, filerule);
    133133        if (!chipTXT) goto finish;
    134134
    … …  
    167167
    168168// extract the input sources corresponding to this readout
    169 // XXX this function needs to be updated to work with the new context of pshot inputs
     169// XXX this function needs to be updated to work with the new context of psphot inputs
    170170psArray *psphotLoadPSFSources (pmConfig *config, const pmFPAview *view) {
    171171
    … …  
    197197// psphotDetectionsFromSources to psphotSourceStats and are now stored on
    198198// detections->newSources.
    199 bool psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view) {
    200 
    201     // find the currently selected readout
    202     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest
     199bool psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view, const char *filerule) {
     200
     201    // find the currently selected readout
     202    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest
    203203    psAssert (file, "missing file?");
    204204
    … …  
    227227// generate the detection structure for the supplied array of sources
    228228// XXX this currently assumes there is a single input file
    229 bool psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, psArray *sources) {
    230 
    231     // find the currently selected readout
    232     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest
     229bool psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *sources) {
     230
     231    // find the currently selected readout
     232    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest
    233233    psAssert (file, "missing file?");
    234234
    … …  
    335335}
    336336
    337 bool psphotCheckExtSources (pmConfig *config, const pmFPAview *view) {
     337bool psphotCheckExtSources (pmConfig *config, const pmFPAview *view, const char *filerule) {
    338338
    339339    bool status;
    … …  
    343343
    344344    // find the currently selected readout
    345     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest
     345    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest
    346346    psAssert (file, "missing file?");
    347347
    … …  
    373373
    374374        // find the detections (by peak and/or footprint) in the image.
    375         if (!psphotFindDetections (config, view, true)) {
     375        if (!psphotFindDetections (config, view, filerule, true)) {
    376376            psError(PSPHOT_ERR_CONFIG, false, "unable to find detections in this image");
    377             return psphotReadoutCleanup (config, view);
     377            return psphotReadoutCleanup (config, view, filerule);
    378378        }
    379379
    380380        // construct sources and measure basic stats
    381         psphotSourceStats (config, view, true);
     381        psphotSourceStats (config, view, filerule, true);
    382382
    383383        // find blended neighbors of very saturated stars
    384         psphotDeblendSatstars (config, view);
     384        psphotDeblendSatstars (config, view, filerule);
    385385
    386386        // mark blended peaks PS_SOURCE_BLEND
    387         if (!psphotBasicDeblend (config, view)) {
     387        if (!psphotBasicDeblend (config, view, filerule)) {
    388388            psLogMsg ("psphot", 3, "failed on deblend analysis");
    389             return psphotReadoutCleanup (config, view);
     389            return psphotReadoutCleanup (config, view, filerule);
    390390        }
    391391
    392392        // classify sources based on moments, brightness
    393         if (!psphotRoughClass (config, view)) {
     393        if (!psphotRoughClass (config, view, filerule)) {
    394394            psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    395             return psphotReadoutCleanup (config, view);
    396         }
    397     }
    398 
    399     return true;
    400 }
     395            return psphotReadoutCleanup (config, view, filerule);
     396        }
     397    }
     398
     399    return true;
     400}
     401
     402// copy the detections from one pmFPAfile to another
     403bool psphotCopySources (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc)
     404{
     405    bool status = true;
     406
     407    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     408    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     409
     410    // skip the chisq image because it is a duplicate of the detection version
     411    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     412    if (!status) chisqNum = -1;
     413
     414    // loop over the available readouts
     415    for (int i = 0; i < num; i++) {
     416        if (i == chisqNum) continue; // skip chisq image
     417        if (!psphotCopySourcesReadout (config, view, ruleOut, ruleSrc, i)) {
     418            psError (PSPHOT_ERR_CONFIG, false, "failed to copy sources from %s to %s entry %d", ruleSrc, ruleOut, i);
     419            return false;
     420        }
     421    }
     422    return true;
     423}
     424
     425// add newly selected sources to the existing list of sources
     426bool psphotCopySourcesReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) {
     427
     428    bool status;
     429
     430    // find the currently selected readout
     431    pmFPAfile *fileSrc = pmFPAfileSelectSingle(config->files, ruleSrc, index); // File of interest
     432    psAssert (fileSrc, "missing file?");
     433
     434    pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa);
     435    psAssert (readoutSrc, "missing readout?");
     436
     437    pmDetections *detections = psMetadataLookupPtr (&status, readoutSrc->analysis, "PSPHOT.DETECTIONS");
     438    psAssert (detections, "missing detections?");
     439
     440    // find the currently selected readout
     441    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, ruleOut, index); // File of interest
     442    psAssert (fileOut, "missing file?");
     443
     444    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
     445    psAssert (readoutOut, "missing readout?");
     446
     447    // save detections on the readout->analysis
     448    // XXX this replaced any existing entry; allow this operation to merge?
     449    if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detections)) {
     450        psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
     451        return false;
     452    }
     453
     454    return true;
     455}
     456

  • trunk/psphot/src/psphotModelBackground.c

    r27657 r28013  
    384384
    385385// XXX supply filename or keep PSPHOT.INPUT fixed?
    386 bool psphotModelBackground (pmConfig *config, const pmFPAview *view)
     386bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule)
    387387{
    388388    bool status = false;
    … …  
    393393    // loop over the available readouts
    394394    for (int i = 0; i < num; i++) {
    395         if (!psphotModelBackgroundReadoutFileIndex(config, view, "PSPHOT.INPUT", i)) {
     395        if (!psphotModelBackgroundReadoutFileIndex(config, view, filerule, i)) {
    396396            psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i);
    397397            return false;

  • trunk/psphot/src/psphotModelTest.c

    r26894 r28013  
    33
    44// XXX add more test information?
    5 bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe) {
     5bool psphotModelTest (pmConfig *config, const pmFPAview *view, const char *filerule, psMetadata *recipe) {
    66
    77    bool status;
    … …  
    3333
    3434    // use poissonian errors or local-sky errors
    35     bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
     35    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, filerule);
    3636    if (!status) POISSON_ERRORS = true;
    3737    pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS);

  • trunk/psphot/src/psphotOutput.c

    r26894 r28013  
    126126}
    127127
    128 bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view) {
     128bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view, const char *filerule) {
    129129
    130130    bool status = false;
    … …  
    135135    // loop over the available readouts
    136136    for (int i = 0; i < num; i++) {
    137         if (!psphotAddPhotcodeReadout (config, view, "PSPHOT.INPUT", i)) {
     137        if (!psphotAddPhotcodeReadout (config, view, filerule, i)) {
    138138            psError (PSPHOT_ERR_CONFIG, false, "failed to add photcode to PSPHOT.INPUT entry %d", i);
    139139            return false;

  • trunk/psphot/src/psphotPetrosianRadialBins.c

    r27819 r28013  
    182182        psFree(values);
    183183        psFree(stats);
    184         return false;
     184        return true;
    185185    }
    186186

  • trunk/psphot/src/psphotPetrosianStats.c

    r27819 r28013  
    1515
    1616    pmSourceRadialProfile *profile = source->extpars->petProfile;
     17
     18    if (!profile->binSB) {
     19        psLogMsg ("psphot", PS_LOG_DETAIL, "no petrosian profile, skipping source %f, %f", source->peak->xf, source->peak->yf);
     20        return true;
     21    }
    1722
    1823    psVector *binSB      = profile->binSB;
    … …  
    190195    source->extpars->petrosianR90Err    = NAN;
    191196
    192     fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
     197    // fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
    193198    psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
    194199

  • trunk/psphot/src/psphotRadialBins.c

    r27819 r28013  
    4444    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    4545    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    46     if (!radMin || !radMin->n) return false;
    47     if (!radMax || !radMax->n) return false;
     46    if (!radMin || !radMin->n) {
     47        psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMin missing or empty)");
     48        return false;
     49    }
     50    if (!radMax || !radMax->n) {
     51        psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMax missing or empty)");
     52        return false;
     53    }
    4854
    4955    psVector *binSB      = psVectorAllocEmpty(radMin->n, PS_TYPE_F32); // surface brightness of radial bin
    … …  
    139145        }
    140146    }
    141     binSB->n = binSBstdev->n = binRad->n = binArea->n = nOut;
     147    binSB->n = binSBstdev->n = binSum->n = binRad->n = binArea->n = nOut;
    142148
    143149    // interpolate any bins that were empty (extrapolate to center if needed)
    … …  
    155161        psFree(values);
    156162        psFree(stats);
    157         return false;
     163        return true;
    158164    }
    159165

  • trunk/psphot/src/psphotReadout.c

    r27657 r28013  
    2828
    2929    // set the photcode for this image
    30     if (!psphotAddPhotcode (config, view)) {
     30    if (!psphotAddPhotcode (config, view, "PSPHOT.INPUT")) {
    3131        psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
    3232        return false;
    … …  
    3434
    3535    // Generate the mask and weight images, including the user-defined analysis region of interest
    36     if (!psphotSetMaskAndVariance (config, view)) {
    37         return psphotReadoutCleanup(config, view);
     36    if (!psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT")) {
     37        return psphotReadoutCleanup(config, view, "PSPHOT.INPUT");
    3838    }
    3939    if (!strcasecmp (breakPt, "NOTHING")) {
    40         return psphotReadoutCleanup(config, view);
     40        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4141    }
    4242
    4343    // generate a background model (median, smoothed image)
    44     if (!psphotModelBackground (config, view)) {
    45         return psphotReadoutCleanup (config, view);
    46     }
    47     if (!psphotSubtractBackground (config, view)) {
    48         return psphotReadoutCleanup (config, view);
     44    if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) {
     45        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
     46    }
     47    if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) {
     48        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4949    }
    5050    if (!strcasecmp (breakPt, "BACKMDL")) {
    51         return psphotReadoutCleanup (config, view);
     51        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    5252    }
    5353
    5454    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are determined and saved on
    5555    // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT
    56     if (!psphotLoadPSF (config, view)) {
     56    if (!psphotLoadPSF (config, view)) { // ??? need to supply 2 ?
    5757        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
    58         return psphotReadoutCleanup (config, view);
     58        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    5959    }
    6060
    6161    // find the detections (by peak and/or footprint) in the image.
    62     if (!psphotFindDetections (config, view, true)) { // pass 1
     62    if (!psphotFindDetections (config, view, "PSPHOT.INPUT", true)) { // pass 1
    6363        // this only happens if we had an error in psphotFindDetections
    6464        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    65         return psphotReadoutCleanup (config, view);
     65        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    6666    }
    6767
    6868    // construct sources and measure basic stats (saved on detections->newSources)
    69     if (!psphotSourceStats (config, view, true)) { // pass 1
     69    if (!psphotSourceStats (config, view, "PSPHOT.INPUT", true)) { // pass 1
    7070        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    71         return psphotReadoutCleanup (config, view);
     71        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    7272    }
    7373    if (!strcasecmp (breakPt, "PEAKS")) {
    74         return psphotReadoutCleanup(config, view);
     74        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    7575    }
    7676
    7777    // find blended neighbors of very saturated stars (detections->newSources)
    78     if (!psphotDeblendSatstars (config, view)) {
     78    if (!psphotDeblendSatstars (config, view, "PSPHOT.INPUT")) {
    7979        psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
    80         return psphotReadoutCleanup (config, view);
     80        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    8181    }
    8282
    8383    // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
    84     if (!psphotBasicDeblend (config, view)) {
     84    if (!psphotBasicDeblend (config, view, "PSPHOT.INPUT")) {
    8585        psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
    86         return psphotReadoutCleanup (config, view);
     86        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    8787    }
    8888
    8989    // classify sources based on moments, brightness.  if a PSF model has been loaded, the PSF
    9090    // clump defined for it is used not measured (detections->newSources)
    91     if (!psphotRoughClass (config, view)) { // pass 1
     91    if (!psphotRoughClass (config, view, "PSPHOT.INPUT")) { // pass 1
    9292        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
    93         return psphotReadoutCleanup (config, view);
     93        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    9494    }
    9595    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
    96     if (!psphotImageQuality (config, view)) { // pass 1
     96    if (!psphotImageQuality (config, view, "PSPHOT.INPUT")) { // pass 1
    9797        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
    98         return psphotReadoutCleanup(config, view);
     98        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    9999    }
    100100    if (!strcasecmp (breakPt, "MOMENTS")) {
    101         return psphotReadoutCleanup(config, view);
     101        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    102102    }
    103103
    104104    // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
    105105    // this step is skipped
    106     if (!psphotChoosePSF (config, view)) { // pass 1
     106    if (!psphotChoosePSF (config, view, "PSPHOT.INPUT")) { // pass 1
    107107        psLogMsg ("psphot", 3, "failure to construct a psf model");
    108         return psphotReadoutCleanup (config, view);
     108        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    109109    }
    110110    if (!strcasecmp (breakPt, "PSFMODEL")) {
    111         return psphotReadoutCleanup (config, view);
     111        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    112112    }
    113113
    114114    // include externally-supplied sources
    115115    // XXX fix this in the new multi-input context
    116     // psphotLoadExtSources (config, view); // pass 1
     116    // psphotLoadExtSources (config, view, "PSPHOT.INPUT"); // pass 1
    117117
    118118    // construct an initial model for each object, set the radius to fitRadius, set circular
    119119    // fit mask (detections->newSources)
    120     psphotGuessModels (config, view); // pass 1
     120    psphotGuessModels (config, view, "PSPHOT.INPUT"); // pass 1
    121121
    122122    // merge the newly selected sources into the existing list
    123123    // NOTE: merge OLD and NEW
    124     psphotMergeSources (config, view);
     124    psphotMergeSources (config, view, "PSPHOT.INPUT");
    125125
    126126    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    127     psphotFitSourcesLinear (config, view, false); // pass 1 (detections->allSources)
     127    psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", false); // pass 1 (detections->allSources)
    128128
    129129    // identify CRs and extended sources (only unmeasured sources are measured)
    130     psphotSourceSize (config, view, true); // pass 1 (detections->allSources)
     130    psphotSourceSize (config, view, "PSPHOT.INPUT", true); // pass 1 (detections->allSources)
    131131    if (!strcasecmp (breakPt, "ENSEMBLE")) {
    132132        goto finish;
    … …  
    135135    // non-linear PSF and EXT fit to brighter sources
    136136    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
    137     psphotBlendFit (config, view); // pass 1 (detections->allSources)
     137    psphotBlendFit (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources)
    138138
    139139    // replace all sources
    140     psphotReplaceAllSources (config, view); // pass 1 (detections->allSources)
     140    psphotReplaceAllSources (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources)
    141141
    142142    // linear fit to include all sources (subtract again)
    143143    // NOTE : apply to ALL sources (extended + psf)
    144     psphotFitSourcesLinear (config, view, true); // pass 2 (detections->allSources)
     144    psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", true); // pass 2 (detections->allSources)
    145145
    146146    // if we only do one pass, skip to extended source analysis
    … …  
    150150
    151151    // add noise for subtracted objects
    152     psphotAddNoise (config, view); // pass 1 (detections->allSources)
     152    psphotAddNoise (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources)
    153153
    154154    // find fainter sources
    155155    // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
    156     psphotFindDetections (config, view, false); // pass 2 (detections->peaks, detections->footprints)
     156    psphotFindDetections (config, view, "PSPHOT.INPUT", false); // pass 2 (detections->peaks, detections->footprints)
    157157
    158158    // remove noise for subtracted objects (ie, return to normal noise level)
    159159    // NOTE: this needs to operate only on the OLD sources
    160     psphotSubNoise (config, view); // pass 1 (detections->allSources)
     160    psphotSubNoise (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources)
    161161
    162162    // define new sources based on only the new peaks
    163163    // NOTE: new sources are saved on detections->newSources
    164     psphotSourceStats (config, view, false); // pass 2 (detections->newSources)
     164    psphotSourceStats (config, view, "PSPHOT.INPUT", false); // pass 2 (detections->newSources)
    165165
    166166    // set source type
    167167    // NOTE: apply only to detections->newSources
    168     if (!psphotRoughClass (config, view)) { // pass 2 (detections->newSources)
     168    if (!psphotRoughClass (config, view, "PSPHOT.INPUT")) { // pass 2 (detections->newSources)
    169169        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    170         return psphotReadoutCleanup (config, view);
     170        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    171171    }
    172172
    173173    // create full input models, set the radius to fitRadius, set circular fit mask
    174174    // NOTE: apply only to detections->newSources
    175     psphotGuessModels (config, view); // pass 2 (detections->newSources)
     175    psphotGuessModels (config, view, "PSPHOT.INPUT"); // pass 2 (detections->newSources)
    176176
    177177    // replace all sources so fit below applies to all at once
    178178    // NOTE: apply only to OLD sources (which have been subtracted)
    179     psphotReplaceAllSources (config, view); // pass 2
     179    psphotReplaceAllSources (config, view, "PSPHOT.INPUT"); // pass 2
    180180
    181181    // merge the newly selected sources into the existing list
    182182    // NOTE: merge OLD and NEW
    183183    // XXX check on free of sources...
    184     psphotMergeSources (config, view); // (detections->newSources + detections->allSources -> detections->allSources)
     184    psphotMergeSources (config, view, "PSPHOT.INPUT"); // (detections->newSources + detections->allSources -> detections->allSources)
    185185
    186186    // NOTE: apply to ALL sources
    187     psphotFitSourcesLinear (config, view, true); // pass 3 (detections->allSources)
     187    psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", true); // pass 3 (detections->allSources)
    188188
    189189pass1finish:
    … …  
    191191    // measure source size for the remaining sources
    192192    // NOTE: applies only to NEW (unmeasured) sources
    193     psphotSourceSize (config, view, false); // pass 2 (detections->allSources)
    194 
    195     psphotExtendedSourceAnalysis (config, view); // pass 1 (detections->allSources)
    196     psphotExtendedSourceFits (config, view); // pass 1 (detections->allSources)
     193    psphotSourceSize (config, view, "PSPHOT.INPUT", false); // pass 2 (detections->allSources)
     194
     195    psphotExtendedSourceAnalysis (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources)
     196    psphotExtendedSourceFits (config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources)
    197197
    198198finish:
    … …  
    202202
    203203    // measure aperture photometry corrections
    204     if (!psphotApResid (config, view)) { // pass 1 (detections->allSources)
     204    if (!psphotApResid (config, view, "PSPHOT.INPUT")) { // pass 1 (detections->allSources)
    205205        psLogMsg ("psphot", 3, "failed on psphotApResid");
    206         return psphotReadoutCleanup (config, view);
     206        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    207207    }
    208208
    209209    // calculate source magnitudes
    210     psphotMagnitudes(config, view); // pass 1 (detections->allSources)
    211 
    212     if (!psphotEfficiency(config, view)) { // pass 1
     210    psphotMagnitudes(config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources)
     211
     212    if (!psphotEfficiency(config, view, "PSPHOT.INPUT")) { // pass 1
    213213        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    214214        psErrorClear();
    … …  
    219219
    220220    // replace background in residual image
    221     psphotSkyReplace (config, view); // pass 1
     221    psphotSkyReplace (config, view, "PSPHOT.INPUT"); // pass 1
    222222
    223223    // drop the references to the image pixels held by each source
    224     psphotSourceFreePixels (config, view); // pass 1
     224    psphotSourceFreePixels (config, view, "PSPHOT.INPUT"); // pass 1
    225225
    226226    // create the exported-metadata and free local data
    227     return psphotReadoutCleanup(config, view);
     227    return psphotReadoutCleanup(config, view, "PSPHOT.INPUT");
    228228}

  • trunk/psphot/src/psphotReadoutCleanup.c

    r27657 r28013  
    22
    33// for now, let's store the detections on the readout->analysis for each readout
    4 bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view)
     4bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view, const char *filerule)
    55{
    66    bool status = true;
    … …  
    2424    // loop over the available readouts
    2525    for (int i = 0; i < num; i++) {
    26         if (!psphotReadoutCleanupReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    27             psError (PSPHOT_ERR_CONFIG, false, "failed on psphotReadoutCleanup for PSPHOT.INPUT entry %d", i);
     26        if (!psphotReadoutCleanupReadout (config, view, filerule, i, recipe)) {
     27            psError (PSPHOT_ERR_CONFIG, false, "failed on psphotReadoutCleanup for %s entry %d", filerule, i);
    2828            return false;
    2929        }
    … …  
    3939// not a DATA error, then there was a serious problem.  Only in this case, or if the fail
    4040// on the stats measurement, do we return false
    41 bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     41bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
    4242
    4343    bool status = true;
    4444
    4545    // find the currently selected readout
    46     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     46    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    4747    psAssert (file, "missing file?");
    4848

  • trunk/psphot/src/psphotReadoutFindPSF.c

    r26894 r28013  
    88
    99    // set the photcode for the PSPHOT.INPUT
    10     if (!psphotAddPhotcode(config, view)) {
     10    if (!psphotAddPhotcode(config, view, "PSPHOT.INPUT")) {
    1111        psError(PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
    1212        return false;
    … …  
    1414
    1515    // Generate the mask and variance images, including the user-defined analysis region of interest
    16     psphotSetMaskAndVariance (config, view);
     16    psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT");
    1717
    1818    // Note that in this implementation, we do NOT model the background and we do not
    … …  
    2121    // include externally-supplied sources (supplied as PSPHOT.INPUT.CMF)
    2222    // XXX we assume a single set of input sources is supplied
    23     if (!psphotDetectionsFromSources (config, view, inSources)) {
     23    if (!psphotDetectionsFromSources (config, view, "PSPHOT.INPUT", inSources)) {
    2424        psError(PSPHOT_ERR_ARGUMENTS, true, "Can't find PSF stars");
    25         return psphotReadoutCleanup(config, view);
     25        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    2626    }
    2727
    2828    // construct detections->newSources and measure basic stats (moments, local sky)
    29     if (!psphotSourceStats(config, view, true)) {
     29    if (!psphotSourceStats(config, view, "PSPHOT.INPUT", true)) {
    3030        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    3131        return false;
    … …  
    3333
    3434    // peak flux is wrong : use the peak measured in the moments analysis:
    35     if (!psphotRepairLoadedSources(config, view)) {
     35    if (!psphotRepairLoadedSources(config, view, "PSPHOT.INPUT")) {
    3636        psError(PSPHOT_ERR_UNKNOWN, false, "failure to repair sources");
    3737        return false;
    … …  
    3939
    4040    // classify sources based on moments, brightness (psf is not known)
    41     if (!psphotRoughClass (config, view)) {
     41    if (!psphotRoughClass (config, view, "PSPHOT.INPUT")) {
    4242        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough source class");
    43         return psphotReadoutCleanup (config, view);
     43        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4444    }
    4545
    46     if (!psphotImageQuality (config, view)) {
     46    if (!psphotImageQuality (config, view, "PSPHOT.INPUT")) {
    4747        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
    48         return psphotReadoutCleanup(config, view);
     48        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4949    }
    5050
    51     if (!psphotChoosePSF(config, view)) {
     51    if (!psphotChoosePSF(config, view, "PSPHOT.INPUT")) {
    5252        psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model");
    53         return psphotReadoutCleanup(config, view);
     53        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    5454    }
    5555
    … …  
    5959    // fits from that analysis, or run the linear PSF fit for all objects currently in hand
    6060    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    61     psphotGuessModels (config, view);
     61    psphotGuessModels (config, view, "PSPHOT.INPUT");
    6262# endif
    6363
    6464    // merge the newly selected sources into the existing list
    6565    // NOTE: merge OLD and NEW
    66     psphotMergeSources (config, view);
     66    psphotMergeSources (config, view, "PSPHOT.INPUT");
    6767
    6868# if 0
    6969    // measure aperture photometry corrections
    70     if (!psphotApResid (config, view)) {
     70    if (!psphotApResid (config, view, "PSPHOT.INPUT")) {
    7171        psLogMsg ("psphot", 3, "failed on psphotApResid");
    72         return psphotReadoutCleanup (config, view);
     72        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    7373    }
    7474# endif
    7575
    7676    // drop the references to the image pixels held by each source
    77     psphotSourceFreePixels(config, view);
     77    psphotSourceFreePixels(config, view, "PSPHOT.INPUT");
    7878
    7979    // create the exported-metadata and free local data
    80     return psphotReadoutCleanup(config, view);
     80    return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    8181}

  • trunk/psphot/src/psphotReadoutKnownSources.c

    r26894 r28013  
    88
    99    // set the photcode for this image
    10     if (!psphotAddPhotcode(config, view)) {
     10    if (!psphotAddPhotcode(config, view, "PSPHOT.INPUT")) {
    1111        psError(PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
    1212        return false;
    … …  
    1414
    1515    // Generate the mask and weight images, including the user-defined analysis region of interest
    16     psphotSetMaskAndVariance (config, view);
     16    psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT");
    1717
    1818    // Note that in this implementation, we do NOT model the background and we do not
    … …  
    2020
    2121    // include externally-supplied sources (supplied as PSPHOT.INPUT.CMF)
    22     if (!psphotDetectionsFromSources (config, view, inSources)) {
     22    if (!psphotDetectionsFromSources (config, view, "PSPHOT.INPUT", inSources)) {
    2323        psError(PSPHOT_ERR_ARGUMENTS, true, "Can't find PSF stars");
    24         return psphotReadoutCleanup(config, view);
     24        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    2525    }
    2626
    2727    // construct sources and measure basic stats
    28     if (!psphotSourceStats (config, view, true)) {
     28    if (!psphotSourceStats (config, view, "PSPHOT.INPUT", true)) {
    2929        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    3030        return false;
    … …  
    3232
    3333    // peak flux is wrong : use the peak measured in the moments analysis:
    34     if (!psphotRepairLoadedSources(config, view)) {
     34    if (!psphotRepairLoadedSources(config, view, "PSPHOT.INPUT")) {
    3535        psError(PSPHOT_ERR_UNKNOWN, false, "failure to repair sources");
    3636        return false;
    … …  
    3838
    3939    // classify sources based on moments, brightness (psf is not known)
    40     if (!psphotRoughClass (config, view)) {
     40    if (!psphotRoughClass (config, view, "PSPHOT.INPUT")) {
    4141        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough source class");
    42         return psphotReadoutCleanup(config, view);
     42        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4343    }
    4444
    45     if (!psphotChoosePSF (config, view)) {
     45    if (!psphotChoosePSF (config, view, "PSPHOT.INPUT")) {
    4646        psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model");
    47         return psphotReadoutCleanup(config, view);
     47        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4848    }
    4949
    5050    // construct an initial model for each object
    51     psphotGuessModels (config, view);
     51    psphotGuessModels (config, view, "PSPHOT.INPUT");
    5252
    5353    // merge the newly selected sources into the existing list
    5454    // NOTE: merge OLD and NEW
    55     psphotMergeSources (config, view);
     55    psphotMergeSources (config, view, "PSPHOT.INPUT");
    5656
    5757    // linear PSF fit to source peaks
    58     psphotFitSourcesLinear (config, view, false);
     58    psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", false);
    5959
    6060    // measure aperture photometry corrections
    61     if (!psphotApResid (config, view)) {
     61    if (!psphotApResid (config, view, "PSPHOT.INPUT")) {
    6262        psLogMsg ("psphot", 3, "failed on psphotApResid");
    63         return psphotReadoutCleanup(config, view);
     63        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    6464    }
    6565
    6666    // calculate source magnitudes
    67     psphotMagnitudes(config, view);
     67    psphotMagnitudes(config, view, "PSPHOT.INPUT");
    6868
    6969    // drop the references to the image pixels held by each source
    70     psphotSourceFreePixels (config, view);
     70    psphotSourceFreePixels (config, view, "PSPHOT.INPUT");
    7171
    7272    // create the exported-metadata and free local data
    73     return psphotReadoutCleanup(config, view);
     73    return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    7474}

  • trunk/psphot/src/psphotReadoutMinimal.c

    r26894 r28013  
    1818
    1919    // set the photcode for this image
    20     if (!psphotAddPhotcode(config, view)) {
     20    if (!psphotAddPhotcode(config, view, "PSPHOT.INPUT")) {
    2121        psError(PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
    2222        return false;
    … …  
    2424
    2525    // Generate the mask and weight images, including the user-defined analysis region of interest
    26     psphotSetMaskAndVariance (config, view);
     26    psphotSetMaskAndVariance (config, view, "PSPHOT.INPUT");
    2727
    2828    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved on readout->analysis
    2929    if (!psphotLoadPSF (config, view)) {
    3030      psError (PSPHOT_ERR_CONFIG, false, "missing psf model");
    31       return psphotReadoutCleanup (config, view);
     31      return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    3232    }
    3333
    3434    // find the detections (by peak and/or footprint) in the image. (final pass)
    35     if (!psphotFindDetections(config, view, false)) {
     35    if (!psphotFindDetections(config, view, "PSPHOT.INPUT", false)) {
    3636        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    37         return psphotReadoutCleanup (config, view);
     37        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    3838    }
    3939
    4040    // construct sources and measure basic stats (saved on detections->newSources)
    41     if (!psphotSourceStats (config, view, false)) { // pass 1
     41    if (!psphotSourceStats (config, view, "PSPHOT.INPUT", false)) { // pass 1
    4242        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    43         return psphotReadoutCleanup (config, view);
     43        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    4444    }
    4545
    4646    // find blended neighbors of very saturated stars
    47     psphotDeblendSatstars (config, view);
     47    psphotDeblendSatstars (config, view, "PSPHOT.INPUT");
    4848
    4949    // mark blended peaks PS_SOURCE_BLEND
    50     if (!psphotBasicDeblend (config, view)) {
     50    if (!psphotBasicDeblend (config, view, "PSPHOT.INPUT")) {
    5151        psLogMsg ("psphot", 3, "failed on deblend analysis");
    52         return psphotReadoutCleanup (config, view);
     52        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    5353    }
    5454
    5555    // classify sources based on moments, brightness (use supplied psf shape parameters)
    56     if (!psphotRoughClass (config, view)) {
     56    if (!psphotRoughClass (config, view, "PSPHOT.INPUT")) {
    5757        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    58         return psphotReadoutCleanup (config, view);
     58        return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    5959    }
    6060
    6161    // construct an initial model for each object
    62     psphotGuessModels (config, view);
     62    psphotGuessModels (config, view, "PSPHOT.INPUT");
    6363
    6464    // merge the newly selected sources into the existing list
    65     psphotMergeSources (config, view);
     65    psphotMergeSources (config, view, "PSPHOT.INPUT");
    6666
    6767    // linear PSF fit to source peaks
    68     psphotFitSourcesLinear (config, view, false);
     68    psphotFitSourcesLinear (config, view, "PSPHOT.INPUT", false);
    6969
    7070// XXX eventually, add the extended source fits here
    7171# if (0)
    7272    // measure source size for the remaining sources
    73     psphotSourceSize (config, view);
     73    psphotSourceSize (config, view, "PSPHOT.INPUT");
    7474
    75     psphotExtendedSourceAnalysis (config, view);
     75    psphotExtendedSourceAnalysis (config, view, "PSPHOT.INPUT");
    7676
    77     psphotExtendedSourceFits (config, view);
     77    psphotExtendedSourceFits (config, view, "PSPHOT.INPUT");
    7878# endif
    7979
    8080    // calculate source magnitudes
    81     psphotMagnitudes(config, view);
     81    psphotMagnitudes(config, view, "PSPHOT.INPUT");
    8282
    8383    // XXX ensure this is measured if the analysis succeeds (even if quality is low)
    84     if (!psphotEfficiency(config, view)) {
     84    if (!psphotEfficiency(config, view, "PSPHOT.INPUT")) {
    8585        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    8686        psErrorClear();
    … …  
    8888
    8989    // drop the references to the image pixels held by each source
    90     psphotSourceFreePixels (config, view);
     90    psphotSourceFreePixels (config, view, "PSPHOT.INPUT");
    9191
    9292    // create the exported-metadata and free local data
    93     return psphotReadoutCleanup(config, view);
     93    return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
    9494}

  • trunk/psphot/src/psphotReplaceUnfit.c

    r26894 r28013  
    2323
    2424// for now, let's store the detections on the readout->analysis for each readout
    25 bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view)
     25bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule)
    2626{
    2727    bool status = true;
    … …  
    3636    // loop over the available readouts
    3737    for (int i = 0; i < num; i++) {
    38         if (!psphotReplaceAllSourcesReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     38        if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe)) {
    3939            psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for PSPHOT.INPUT entry %d", i);
    4040            return false;

  • trunk/psphot/src/psphotRoughClass.c

    r27657 r28013  
    88
    99// for now, let's store the detections on the readout->analysis for each readout
    10 bool psphotRoughClass (pmConfig *config, const pmFPAview *view)
     10bool psphotRoughClass (pmConfig *config, const pmFPAview *view, const char *filerule)
    1111{
    1212    bool status = true;
    … …  
    2626    for (int i = 0; i < num; i++) {
    2727        if (i == chisqNum) continue; // skip chisq image
    28         if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    29             psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for PSPHOT.INPUT entry %d", i);
     28        if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) {
     29            psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for %s entry %d", filerule, i);
    3030            return false;
    3131        }
    … …  
    3434}
    3535
    36 bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     36bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
    3737
    3838    bool status;
    … …  
    4141
    4242    // find the currently selected readout
    43     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     43    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    4444    psAssert (file, "missing file?");
    4545

  • trunk/psphot/src/psphotSetMaskBits.c

    r21254 r28013  
    3737    return true;
    3838}
     39
     40// XXX should these be in config->analysis or somewhere else besides 'recipe'?

  • trunk/psphot/src/psphotSkyReplace.c

    r27657 r28013  
    11# include "psphotInternal.h"
    22
    3 bool psphotSkyReplace (pmConfig *config, const pmFPAview *view)
     3bool psphotSkyReplace (pmConfig *config, const pmFPAview *view, const char *filerule)
    44{
    55    bool status = true;
    … …  
    1515    for (int i = 0; i < num; i++) {
    1616        if (i == chisqNum) continue; // skip chisq image
    17         if (!psphotSkyReplaceReadout (config, view, "PSPHOT.INPUT", i)) {
    18             psError (PSPHOT_ERR_CONFIG, false, "failed to replace sky for PSPHOT.INPUT entry %d", i);
     17        if (!psphotSkyReplaceReadout (config, view, filerule, i)) {
     18            psError (PSPHOT_ERR_CONFIG, false, "failed to replace sky for %s entry %d", filerule, i);
    1919            return false;
    2020        }
    … …  
    2525// XXX make this an option?
    2626// in order to  successfully replace the sky, we must define a corresponding file...
    27 bool psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
     27bool psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
    2828
    2929    psTimerStart ("psphot.skyreplace");
    3030
    3131    // find the currently selected readout
    32     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     32    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    3333    psAssert (file, "missing file?");
    3434

  • trunk/psphot/src/psphotSourceFreePixels.c

    r26894 r28013  
    11# include "psphotInternal.h"
    22
    3 bool psphotSourceFreePixels (pmConfig *config, const pmFPAview *view)
     3bool psphotSourceFreePixels (pmConfig *config, const pmFPAview *view, const char *filerule)
    44{
    55    bool status = true;
    … …  
    1010    // loop over the available readouts
    1111    for (int i = 0; i < num; i++) {
    12         if (!psphotSourceFreePixelsReadout (config, view, "PSPHOT.INPUT", i)) {
    13             psError (PSPHOT_ERR_CONFIG, false, "failed to free source pixels for PSPHOT.INPUT entry %d", i);
     12        if (!psphotSourceFreePixelsReadout (config, view, filerule, i)) {
     13            psError (PSPHOT_ERR_CONFIG, false, "failed to free source pixels for %s entry %d", filerule, i);
    1414            return false;
    1515        }
    … …  
    1818}
    1919
    20 bool psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index) {
     20bool psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
    2121
    2222    bool status;
    2323
    2424    // find the currently selected readout
    25     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     25    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    2626    psAssert (file, "missing file?");
    2727

  • trunk/psphot/src/psphotSourceMatch.c

    r27657 r28013  
    11# include "psphotInternal.h"
    22
    3 bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects);
    4  
    5 psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view)
     3bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
     4bool psphotMatchSourcesSetIDs (psArray *objects);
     5 
     6psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view, const char *filerule)
    67{
    78    bool status = true;
    … …  
    1415    // loop over the available readouts
    1516    for (int i = 0; i < num; i++) {
    16         if (!psphotMatchSourcesReadout (objects, config, view, "PSPHOT.INPUT", i)) {
     17        if (!psphotMatchSourcesReadout (objects, config, view, filerule, i)) {
    1718            psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i);
    1819            psFree (objects);
    … …  
    2122    }
    2223
    23     psphotMatchSourcesGenerate (config, view, objects);
     24    // create sources for images where an object has been detected in the other images
     25    psphotMatchSourcesAddMissing (config, view, filerule, objects);
     26
     27    // choose a consistent position; set common sequence values
     28    psphotMatchSourcesSetIDs (objects);
    2429
    2530    return objects;
    2631}
    2732
    28 bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index) {
     33bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
    2934 
    3035    bool status = false;
    … …  
    3843
    3944    // find the currently selected readout
    40     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     45    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    4146    psAssert (file, "missing file?");
    4247
    … …  
    145150}
    146151
    147 bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects) {
     152bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
    148153
    149154    bool status = false;
    … …  
    167172
    168173        // find the currently selected readout
    169         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest
     174        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
    170175        psAssert (file, "missing file?");
    171176
    … …  
    255260    return true;
    256261}
     262
     263bool psphotMatchSourcesSetIDs (psArray *objects) {
     264
     265    for (int i = 0; i < objects->n; i++) {
     266        pmPhotObj *obj = objects->data[i];
     267
     268        // set the source->seq values
     269        for (int j = 0; j < obj->sources->n; j++) {
     270            pmSource *src = obj->sources->data[j];
     271            src->seq = i;
     272        }
     273    }
     274    return true;
     275}

  • trunk/psphot/src/psphotSourceSize.c

    r27695 r28013  
    3333
    3434// for now, let's store the detections on the readout->analysis for each readout
    35 bool psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize)
     35bool psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize)
    3636{
    3737    bool status = true;
    … …  
    5151    for (int i = 0; i < num; i++) {
    5252        if (i == chisqNum) continue; // skip chisq image
    53         if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) {
     53        if (!psphotSourceSizeReadout (config, view, filerule, i, recipe, getPSFsize)) {
    5454            psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i);
    5555            return false;
    … …  
    6060
    6161// this function use an internal flag to mark sources which have already been measured
    62 bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool getPSFsize)
     62bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool getPSFsize)
    6363{
    6464    bool status;
    … …  
    6868
    6969    // find the currently selected readout
    70     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     70    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    7171    psAssert (file, "missing file?");
    7272

  • trunk/psphot/src/psphotSourceStats.c

    r27657 r28013  
    55// The new sources are added to any existing sources on detections->newSources.  The sources
    66// on detections->allSources are ignored.
    7 bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow)
     7bool psphotSourceStats (pmConfig *config, const pmFPAview *view, const char *filerule, bool setWindow)
    88{
    99    bool status = true;
    … …  
    1616    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    1717
     18    // skip the chisq image (optionally?)
     19    // int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     20    // if (!status) chisqNum = -1;
     21
    1822    // loop over the available readouts
    1923    for (int i = 0; i < num; i++) {
    20         if (!psphotSourceStatsReadout (config, view, "PSPHOT.INPUT", i, recipe, setWindow)) {
    21             psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
     24        // if (i == chisqNum) continue; // skip chisq image
     25        if (!psphotSourceStatsReadout (config, view, filerule, i, recipe, setWindow)) {
     26            psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for %s entry %d", filerule, i);
    2227            return false;
    2328        }
    … …  
    2631}
    2732
    28 bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool setWindow) {
     33bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool setWindow) {
    2934
    3035    bool status = false;
    … …  
    3439
    3540    // find the currently selected readout
    36     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     41    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    3742    psAssert (file, "missing file?");
    3843

  • trunk/psphot/src/psphotStackArguments.c

    r27657 r28013  
    11# include "psphotStandAlone.h"
    22
     3static void dumpTemplate(void);
    34static void usage(pmConfig *config, int exitCode);
    45static void writeHelpInfo(FILE* ofile);
    … …  
    1112    if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(stdout);
    1213    if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(stdout);
     14
     15    if (psArgumentGet(argc, argv, "-template")) dumpTemplate();
    1316
    1417    // load config data from default locations
    … …  
    8487          "\n"
    8588          "where INPUTS.mdc contains various METADATAs, each with:\n"
    86           "\tIMAGE(STR):     Image filename\n"
    87           "\tMASK(STR):      Mask filename\n"
    88           "\tVARIANCE(STR):  Variance map filename\n"
     89          "\tIMAGE    : Image filename\n"
     90          "\tMASK     : Mask filename\n"
     91          "\tVARIANCE : Variance map filename\n"
     92          "(use -template to generate a sample input.mdc file)\n"
    8993          "OUTROOT is the 'root name' for output files\n"
    9094          "\n"
    … …  
    144148}
    145149
     150static void dumpTemplate(void) {
     151
     152    fprintf (stdout, "# this line is required for multiple INPUT blocks to be accepted\n");
     153    fprintf (stdout, "INPUT MULTI\n\n");
     154
     155    fprintf (stdout, "# copy and repeat the following block as needed (one per input image set)\n");
     156    fprintf (stdout, "# either RAW (unconvolved) or CNV (convolved) input images are required\n");
     157    fprintf (stdout, "# if both are supplied, by default RAW is used for detection, CNV is convolved (further) to match target PSF\n");
     158    fprintf (stdout, "# if MASK or VARIANCE images are not supplied, they will be generated\n");
     159    fprintf (stdout, "# if MASK or VARIANCE images are not supplied, they will be generated\n");
     160    fprintf (stdout, "# PSF may be supplied for the convolution target\n");
     161    fprintf (stdout, "INPUT METADATA\n");
     162    fprintf (stdout, "  RAW:IMAGE     STR   file.im.fits   # signal image filename\n");
     163    fprintf (stdout, "  RAW:MASK      STR   file.mk.fits   # mask image filename\n");
     164    fprintf (stdout, "  RAW:VARIANCE  STR   file.wt.fits   # variance image filename\n");
     165    fprintf (stdout, "  RAW:PSF       STR   file.psf.fits  # psf from input unconvolved image\n");
     166
     167    fprintf (stdout, "  CNV:IMAGE     STR   file.im.fits   # signal image filename\n");
     168    fprintf (stdout, "  CNV:MASK      STR   file.mk.fits   # mask image filename\n");
     169    fprintf (stdout, "  CNV:VARIANCE  STR   file.wt.fits   # variance image filename\n");
     170    fprintf (stdout, "  CNV:PSF       STR   file.psf.fits  # psf from input convolved image\n");
     171
     172    fprintf (stdout, "  SOURCES       STR   file.cmf       # measured source positions\n");
     173    fprintf (stdout, "END\n");
     174
     175    psLibFinalize();
     176    exit(PS_EXIT_SUCCESS);
     177}

  • trunk/psphot/src/psphotStackChisqImage.c

    r27657 r28013  
    66
    77// XXX supply filename or keep PSPHOT.INPUT fixed?
    8 bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view)
     8bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleCnv)
    99{
    1010    bool status = false;
    … …  
    2121
    2222    // loop over the available readouts
     23    // generate the chisq image from the 'detection' images
    2324    for (int i = 0; i < num; i++) {
    24         if (!psphotStackChisqImageAddReadout(config, view, &chiReadout, "PSPHOT.INPUT", i)) {
    25             psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i);
     25        if (!psphotStackChisqImageAddReadout(config, view, ruleDet, &chiReadout, i)) {
     26            psError (PSPHOT_ERR_CONFIG, false, "failed to model background for %s entry %d", ruleDet, i);
    2627            return false;
    2728        }
    … …  
    3536    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", num);
    3637
    37     // need to save the resulting image somewhere (PSPHOT.INPUT?)
    38     if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
     38    // save the resulting image in the 'detection' set
     39    if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleDet, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
     40        psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files");
     41        return false;
     42    }
     43
     44    // save the resulting image in the 'convolved' set
     45    if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleCnv, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
    3946        psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files");
    4047        return false;
    … …  
    4855bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration
    4956                                     const pmFPAview *view,
     57                                     const char *filerule,
    5058                                     pmReadout **chiReadout,
    51                                      char *filename,
    5259                                     int index)
    5360{
    … …  
    5562
    5663    // find the currently selected readout
    57     pmFPAfile *input = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     64    pmFPAfile *input = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    5865    psAssert (input, "missing file?");
    5966
    … …  
    111118}
    112119
    113 bool psphotStackRemoveChisqFromInputs (pmConfig *config) {
     120bool psphotStackRemoveChisqFromInputs (pmConfig *config, const char *filerule) {
    114121
    115122    bool status = false;
    … …  
    121128    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    122129
    123     pmFPAfileRemoveSingle (config->files, "PSPHOT.INPUT", chisqNum);
     130    pmFPAfileRemoveSingle (config->files, filerule, chisqNum);
    124131
    125132    inputNum --;

  • trunk/psphot/src/psphotStackImageLoop.c

    r27848 r28013  
    11# include "psphotStandAlone.h"
    2 
    3 # define ESCAPE(MESSAGE) { \
    4   psError(PSPHOT_ERR_DATA, false, MESSAGE); \
    5   psFree (view); \
    6   return false; \
    7 }
     2#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
     3
     4# define ESCAPE(MESSAGE) {                              \
     5        psError(PSPHOT_ERR_DATA, false, MESSAGE);       \
     6        psFree (view);                                  \
     7        return false;                                   \
     8    }
     9
     10bool GetAstrometryFPA (pmConfig *config, pmFPAview *view);
     11bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry);
     12bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry);
    813
    914bool psphotStackImageLoop (pmConfig *config) {
    … …  
    1419    pmReadout *readout;
    1520
    16     pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
    17     if (!status) {
     21    pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
     22    pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
     23    pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
     24
     25    if (!input) {
    1826        psError(PSPHOT_ERR_PROG, false, "Can't find input data!");
    1927        return false;
    … …  
    2432    // XXX for now, just load the full set of images up front
    2533    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
     34
     35    bool bilevelAstrometry = GetAstrometryFPA(config, view);
    2636
    2737    // for psphot, we force data to be read at the chip level
    … …  
    4151                if (! readout->data_exists) { continue; }
    4252
    43 # if (0)               
    44                 // uncomment to generate matched psfs
     53                // PSF matching
    4554                if (!psphotStackMatchPSFs (config, view)) {
    4655                    psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
    … …  
    4857                    return false;
    4958                }
    50 # endif
    5159
    5260                // XXX for now, we assume there is only a single chip in the PHU:
    … …  
    6977            }
    7078        }
     79
     80        GetAstrometryChip(config, view, bilevelAstrometry);
     81
    7182        // save output which is saved at the chip level
    7283        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot.");
    7384    }
     85
     86    SetAstrometryFPA(config, view, bilevelAstrometry);
     87   
    7488    // save output which is saved at the fpa level
    7589    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot.");
    … …  
    87101*/
    88102
     103# define UPDATE_HEADER 0
     104
     105bool GetAstrometryFPA (pmConfig *config, pmFPAview *view) {
     106
     107    bool status = false;
     108
     109    bool foundAstrometry = false;
     110    bool bilevelAstrometry = false;
     111
     112    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     113    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     114
     115    // loop over the available readouts
     116    for (int i = 0; i < num; i++) {
     117
     118        // find the currently selected readout
     119        pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest
     120        psAssert (output, "missing file?");
     121
     122        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
     123        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
     124        pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
     125        psAssert (input, "missing input file");
     126
     127        // find the FPA phu
     128        pmHDU *phu = pmFPAviewThisPHU(view, input->fpa);
     129        if (!phu) {
     130            psWarning("Unable to read bilevel mosaic astrometry for input FPA entry %d", i);
     131            continue;
     132        }
     133
     134        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
     135        if (!ctype) {
     136            psWarning("Error in WCS keywords for input FPA entry %d", i);
     137            continue;
     138        }
     139
     140        if (!foundAstrometry) {
     141            bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
     142        } else {
     143            if (bilevelAstrometry != !strcmp (&ctype[4], "-DIS")) {
     144                psAbort("astrometry type mis-match");
     145            }
     146        }
     147
     148        if (bilevelAstrometry) {
     149            // update the output structures
     150            if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) {
     151                psWarning("Unable to read bilevel mosaic astrometry for input FPA.");
     152            }
     153        }
     154    }
     155    return bilevelAstrometry;
     156}
     157
     158bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) {
     159
     160    bool status = false;
     161
     162    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     163    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     164
     165    // loop over the available readouts
     166    for (int i = 0; i < num; i++) {
     167
     168        // find the currently selected readout
     169        pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest
     170        psAssert (output, "missing file?");
     171
     172        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
     173        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
     174        pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
     175        psAssert (input, "missing input file");
     176
     177        // Need to update the output for astrometry.  Read WCS data from the corresponding
     178        // header and save in the output fpa
     179        pmHDU *inHDU = pmFPAviewThisHDU (view, input->fpa);
     180        pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
     181
     182# if (UPDATE_HEADER)
     183        pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa);
     184        if (!outHDU) {
     185            pmFPAAddSourceFromView(output->fpa, "name", view, output->format);
     186            outHDU = pmFPAviewThisHDU (view, output->fpa);
     187            psAssert (outHDU, "failed to make HDU");
     188        }
     189# endif
     190
     191        if (bilevelAstrometry) {
     192            if (!pmAstromReadBilevelChip (outChip, inHDU->header)) {
     193                psWarning("Unable to read bilevel chip astrometry for input FPA.");
     194                continue;
     195            }
     196# if (UPDATE_HEADER)
     197            if (!pmAstromWriteBilevelChip(outHDU->header, outChip, WCS_NONLIN_TOL)) {
     198                psWarning("Unable to generate WCS header.");
     199                continue;
     200            }
     201# endif
     202        } else {
     203            // we use a default FPA pixel scale of 1.0
     204            if (!pmAstromReadWCS (output->fpa, outChip, inHDU->header, 1.0)) {
     205                psWarning("Unable to read WCS astrometry for input FPA.");
     206                continue;
     207            }
     208# if (UPDATE_HEADER)
     209            if (UPDATE_HEADER && !pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) {
     210                psWarning("Unable to generate WCS header.");
     211                continue;
     212            }
     213# endif
     214        }
     215    }
     216
     217    return true;
     218}
     219
     220bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) {
     221
     222    bool status = false;
     223
     224    if (!bilevelAstrometry) return true;
     225
     226    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     227    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     228
     229    // loop over the available readouts
     230    for (int i = 0; i < num; i++) {
     231
     232        // find the currently selected readout
     233        pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest
     234        psAssert (output, "missing file?");
     235
     236# if (UPDATE_HEADER)
     237        pmHDU *PHU = pmFPAviewThisPHU(view, output->fpa);
     238        if (!PHU) {
     239            pmFPAAddSourceFromView(output->fpa, "name", view, output->format);
     240            PHU = pmFPAviewThisPHU (view, output->fpa);
     241            psAssert (PHU, "failed to make PHU");
     242        }
     243
     244        if (!pmAstromWriteBilevelMosaic(PHU->header, output->fpa, WCS_NONLIN_TOL)) {
     245            psWarning("Unable to generate WCS header.");
     246        }
     247# endif
     248    }
     249
     250    return true;
     251}

  • trunk/psphot/src/psphotStackMatchPSFs.c

    r27850 r28013  
    11# include "psphotInternal.h"
    22
    3 bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view, bool firstPass)
     3bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view)
    44{
    55    bool status = true;
     6
     7    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PSPHOT");
     8    psAssert(recipe, "We've thrown an error on this before.");
    69
    710    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    811    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    912
    10     // loop over the available readouts
     13    // 'options' carries info needed to perform the stack matching
     14    psphotStackOptions *options = psphotStackOptionsAlloc(num);
     15
     16    options->convolve = psMetadataLookupBool (&status, recipe, "PSPHOT.STACK.MATCH.PSF");
     17    psAssert (status, "PSPHOT.STACK.MATCH.PSF not in recipe");
     18
     19    if (options->convolve) {
     20        char *convolveSource = psMetadataLookupStr (&status, recipe, "PSPHOT.STACK.MATCH.PSF.SOURCE");
     21        options->convolveSource = psphotStackConvolveSourceFromString (convolveSource);
     22        if (options->convolveSource == PSPHOT_CNV_SRC_NONE) {
     23            psError (PSPHOT_ERR_CONFIG, true, "stack convolution source not defined in recipe");
     24            return false;
     25        }
     26    }
     27
     28    // loop over the available readouts (ignore chisq image)
    1129    for (int i = 0; i < num; i++) {
    12         if (!psphotStackMatchPSFsReadout (config, view, i)) {
     30        if (!psphotStackMatchPSFsPrepare (config, view, options, i)) {
     31            psError (PSPHOT_ERR_CONFIG, false, "failed to set PSF matching options for entry %d", i);
     32            return false;
     33        }
     34    }
     35
     36    // Generate target PSF
     37    if (options->convolve) {
     38        options->psf = psphotStackPSF(config, options->numCols, options->numRows, options->psfs, options->inputMask);
     39        if (!options->psf) {
     40            psError(psErrorCodeLast(), false, "Unable to determine output PSF.");
     41            return false;
     42        }
     43        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
     44        options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * options->numCols, 0.5 * options->numRows); // FWHM for target
     45        psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);
     46
     47        // XXX is this needed to supply the psf to psphot??
     48        // pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip
     49        // psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "Target PSF", options->psf);
     50        // outChip->data_exists = true;
     51    }
     52
     53    // loop over the available readouts (ignore chisq image)
     54    for (int i = 0; i < num; i++) {
     55        if (!psphotStackMatchPSFsReadout (config, view, options, i)) {
    1356            psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
    1457            return false;
    1558        }
    1659    }
     60
     61    psFree (options);
    1762    return true;
    1863}
    1964
    2065// convolve the image to match desired PSF
    21 bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, int index) {
     66bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) {
    2267
    23     bool status;
    24     int pass;
    25     float NSIGMA_PEAK = 25.0;
    26     int NMAX = 0;
     68    pmFPAfile *fileSrc = psphotStackGetConvolveSource(config, options, index);
     69    if (!fileSrc) {
     70        psError(PSPHOT_ERR_CONFIG, false, "desired convolution source is missing");
     71        return false;
     72    }
    2773
    28     // find the currently selected readout
    29     pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
    30     psAssert (fileRaw, "missing file?");
     74    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", index); // File of interest
     75    psAssert (fileOut, "missing output file?");
    3176
    32     pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT.CONV", index); // File of interest
    33     psAssert (fileCnv, "missing file?");
     77    pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa);
     78    psAssert (readoutSrc, "missing readout?");
    3479
    35     pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa);
    36     psAssert (readoutRaw, "missing readout?");
     80    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
     81    if (readoutOut == NULL) {
     82        readoutOut = pmFPAGenerateReadout(config, view, "PSPHOT.STACK.OUTPUT.IMAGE", fileSrc->fpa, NULL, index);
     83        psAssert (readoutOut, "missing readout?");
     84    }
    3785
    38     pmReadout *readoutCnv = pmFPAviewThisReadout(view, fileCnv->fpa);
    39     psAssert (readoutCnv, "missing readout?");
    40 
    41     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    42     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    43     psAssert (maskVal, "missing mask value?");
    44 
    45     /***** set up recipe options *****/
    46 
    47     psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    48     psAssert(stackRecipe, "We've thrown an error on this before.");
    49 
    50     // Look up appropriate values from the ppSub recipe
    51     psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    52     psAssert(subRecipe, "recipe missing");
    53 
    54     int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size
    55 
    56     float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
    57 
    58     psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
    59     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    60     psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor
    61     psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
    62     psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
    63     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    64 
    65     bool mdok;                          // Status of MD lookup
    66     float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness
    67     int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
    68 
    69     if (!pmReadoutMaskNonfinite(readout, maskVal)) {
     86    // set NAN pixels to 'SAT'
     87    psImageMaskType maskVal = pmConfigMaskGet("SAT", config);
     88    if (!pmReadoutMaskNonfinite(readoutSrc, maskVal)) {
    7089        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout.");
    7190        return false;
    … …  
    7493    // Image Matching (PSFs or just flux)
    7594    if (options->convolve) {
    76       // Full match of PSFs
    77         pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily
    78 
    79         // Read previously produced kernel
    80         if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
    81           loadKernel();
    82         } else {
    83           matchKernel();
    84         } // !DEBUG.STACK
    85 
    86         saveMatchData();
    87 
    88         saveChiSquare();
    89 
    90         renormKernel();
    91 
    92         // Reject image completely if the maximum deconvolution fraction exceeds the limit
    93         float deconv = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
    94         if (deconv > deconvLimit) {
    95             psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index);
    96             psFree(conv);
    97             return NULL;
    98         }
    99 
    100         readout->analysis = psMetadataCopy(readout->analysis, conv->analysis);
    101 
    102         psFree(conv);
     95        matchKernel(config, readoutOut, readoutSrc, options, index);
     96        saveMatchData(readoutOut, options, index);
     97        // renormKernel(readoutCnv, options, index);
    10398    } else {
    104         // only match the flux
    105         float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
    106         psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
    107         psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
     99        // only match the flux (NO! not for multi-filter, at least!)
     100        // XXX do not generate readoutCnv in this case?
     101        // float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
     102        // psBinaryOp(readoutRaw->image, readoutRaw->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
     103        // psBinaryOp(readoutRaw->variance, readoutRaw->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
    108104    }
    109105
    110     // Ensure the background value is zero
    111     psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
    112     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    113     if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
    114         psWarning("Can't measure background for image.");
    115         psErrorClear();
    116     } else {
    117         if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
    118             psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
    119                      psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
    120             (void)psBinaryOp(readout->image, readout->image, "-",
    121                              psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
    122         }
    123     }
     106    rescaleData(readoutOut, config, options, index);
    124107
    125     if (!stackRenormaliseReadout(config, readout)) {
    126         psFree(rng);
    127         psFree(bg);
    128         return false;
    129     }
    130 
    131     // Measure the variance level for the weighting
    132     if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) {
    133         if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
    134             psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image.");
    135             psFree(rng);
    136             psFree(bg);
    137             return false;
    138         }
    139         options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance));
    140     } else {
    141         options->weightings->data.F32[index] = 1.0;
    142     }
    143     psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n",
    144              index, options->weightings->data.F32[index]);
    145 
    146     psFree(rng);
    147     psFree(bg);
    148 
    149     dumpImage3();
     108    dumpImage(readoutOut, readoutSrc, index, "convolved");
    150109
    151110    return true;
    152111}
     112
     113
     114# if (0)
     115// Read previously produced kernel
     116if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
     117    loadKernel(config, readoutCnv, options, index);
     118} else {
     119    matchKernel(config, readoutCnv, readoutRaw, options, index);
     120}
     121# endif

  • trunk/psphot/src/psphotStackMatchPSFsUtils.c

    r27850 r28013  
    1 /***** defines *****/
    2 
    3 #define ARRAY_BUFFER 16                 // Number to add to array at a time
    4 #define MAG_IGNORE 50                   // Ignore magnitudes fainter than this --- they're not real!
    5 #define FAKE_SIZE 1                     // Size of fake convolution kernel
    6 #define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
    7 #define NOISE_FRACTION 0.01             // Set minimum flux to this fraction of noise
    8 #define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
     1# include "psphotInternal.h"
     2# define ARRAY_BUFFER 16                 // Number to add to array at a time
    93
    104// XXX better name
    … …  
    1812    psFree(resolved);
    1913    if (!fits) {
    20         psError(PPSTACK_ERR_IO, false, "Unable to open previously produced image: %s", name);
     14        psError(PSPHOT_ERR_IO, false, "Unable to open previously produced image: %s", name);
    2115        return false;
    2216    }
    2317    psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
    2418    if (!image) {
    25         psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image: %s", name);
     19        psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image: %s", name);
    2620        psFitsClose(fits);
    2721        return false;
    … …  
    5145}
    5246
     47# define SN_MIN 50.0
    5348psArray *stackSourcesFilter(psArray *sources, // Source list to filter
    54                                    int exclusion // Exclusion zone, pixels
     49                            int exclusion // Exclusion zone, pixels
    5550    )
    5651{
    … …  
    6863            continue;
    6964        }
     65        if (!source->peak) continue;
     66        if (source->peak->SN < SN_MIN) continue;
    7067        coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source);
    7168        numGood++;
    … …  
    8380            continue;
    8481        }
     82        if (!source->peak) continue;
     83        if (source->peak->SN < SN_MIN) continue;
    8584        float xSource, ySource;         // Coordinates of source
    8685        coordsFromSource(&xSource, &ySource, source);
    … …  
    9089
    9190        long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
    92         psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
     91        psTrace("psphotStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
    9392                coords->data.F64[0], coords->data.F64[1], numWithin);
    9493        if (numWithin == 1) {
    … …  
    104103    psFree(y);
    105104
    106     psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
     105    psLogMsg("psphotStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
    107106
    108107    return filtered;
    … …  
    111110// Add background into the fake image
    112111// Based on ppSubBackground()
    113 static psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model
     112psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model
    114113                                     const pmConfig *config // Configuration
    115114    )
    … …  
    121120    int numCols = image->numCols, numRows = image->numRows; // Size of image
    122121
    123     psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE);
     122    psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, "PPSTACK");
    124123    psAssert(ppStackRecipe, "Need PPSTACK recipe");
    125     psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE);
     124    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, "PSPHOT");
    126125    psAssert(psphotRecipe, "Need PSPHOT recipe");
    127126
    … …  
    138137    psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model
    139138    if (!psImageUnbin(unbinned, binned, binning)) {
    140         psError(PPSTACK_ERR_DATA, false, "Unable to unbin background model");
     139        psError(PSPHOT_ERR_DATA, false, "Unable to unbin background model");
    141140        psFree(binned);
    142141        psFree(unbinned);
    … …  
    153152    )
    154153{
    155 #if 1
    156154    bool mdok; // Status of metadata lookups
    157155
    158     psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppStack
     156    psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, "PPSTACK"); // Recipe for ppStack
    159157    psAssert(recipe, "Need PPSTACK recipe");
    160158
    … …  
    163161    int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM");
    164162    if (!mdok) {
    165         psError(PPSTACK_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe");
     163        psError(PSPHOT_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe");
    166164        return false;
    167165    }
    168166    float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN");
    169167    if (!mdok) {
    170         psError(PPSTACK_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe");
     168        psError(PSPHOT_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe");
    171169        return false;
    172170    }
    173171    float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX");
    174172    if (!mdok) {
    175         psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe");
     173        psError(PSPHOT_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe");
    176174        return false;
    177175    }
    … …  
    181179    psImageCovarianceTransfer(readout->variance, readout->covariance);
    182180    return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
    183 #else
    184     return true;
    185 #endif
    186181}
    187182
    … …  
    190185// It implicitly assumes the output root name is the same between invocations.
    191186
    192 bool loadKernel () {
    193             pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
    194             psAssert(file, "Require file");
    195 
    196             pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
    197             view->chip = view->cell = view->readout = 0;
    198             psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest
    199 
    200             // Read convolution kernel
    201             psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
    202             psFree(filename);
    203             psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
    204             psFree(resolved);
    205             if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) {
    206                 psError(PPSTACK_ERR_IO, false, "Unable to read previously produced kernel");
    207                 psFitsClose(fits);
    208                 return false;
    209             }
    210             psFitsClose(fits);
    211 
    212             if (!readImage(&readout->image, options->convImages->data[index], config) ||
    213                 !readImage(&readout->mask, options->convMasks->data[index], config) ||
    214                 !readImage(&readout->variance, options->convVariances->data[index], config)) {
    215                 psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image.");
    216                 return false;
    217             }
    218 
    219             psRegion *region = psMetadataLookupPtr(NULL, conv->analysis,
    220                                                    PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
    221             pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis,
    222                                                                 PM_SUBTRACTION_ANALYSIS_KERNEL);
    223 
    224             pmSubtractionAnalysis(conv->analysis, NULL, kernels, region,
    225                                   readout->image->numCols, readout->image->numRows);
    226 
    227             psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
    228             bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
    229             psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
    230             psImageCovarianceSetThreads(oldThreads);
    231             psFree(readout->covariance);
    232             readout->covariance = covar;
    233             psFree(kernel);
    234 }
    235 
    236 bool dumpImage() {
    237     // XXX should be optional
    238             {
    239                 pmHDU *hdu = pmHDUFromCell(readout->parent);
    240                 psString name = NULL;
    241                 psStringAppend(&name, "fake_%03d.fits", index);
    242                 pmStackVisualPlotTestImage(fake->image, name);
    243                 psFits *fits = psFitsOpen(name, "w");
    244                 psFree(name);
    245                 psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
    246                 psFitsClose(fits);
    247             }
    248             {
    249                 pmHDU *hdu = pmHDUFromCell(readout->parent);
    250                 psString name = NULL;
    251                 psStringAppend(&name, "real_%03d.fits", index);
    252                 pmStackVisualPlotTestImage(readout->image, name);
    253                 psFits *fits = psFitsOpen(name, "w");
    254                 psFree(name);
    255                 psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
    256                 psFitsClose(fits);
    257             }
    258 }
    259 
    260 bool dumpImage2() {
    261     // XXX should be optional
    262 
    263             {
    264                 pmHDU *hdu = pmHDUFromCell(readout->parent);
    265                 psString name = NULL;
    266                 psStringAppend(&name, "conv_%03d.fits", index);
    267                 pmStackVisualPlotTestImage(conv->image, name);
    268                 psFits *fits = psFitsOpen(name, "w");
    269                 psFree(name);
    270                 psFitsWriteImage(fits, hdu->header, conv->image, 0, NULL);
    271                 psFitsClose(fits);
    272             }
    273             {
    274                 pmHDU *hdu = pmHDUFromCell(readout->parent);
    275                 psString name = NULL;
    276                 psStringAppend(&name, "diff_%03d.fits", index);
    277                 pmStackVisualPlotTestImage(fake->image, name);
    278                 psFits *fits = psFitsOpen(name, "w");
    279                 psFree(name);
    280                 psBinaryOp(fake->image, conv->image, "-", fake->image);
    281                 psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
    282                 psFitsClose(fits);
    283             }
    284 }
    285 
    286 bool dumpImage3()
     187# if (0)
     188bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index) {
     189
     190    // Read the convolution kernel from the saved file
     191    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
     192    psAssert(file, "Require file");
     193
     194    pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
     195    view->chip = view->cell = view->readout = 0;
     196    psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest
     197
     198    // Read convolution kernel data
     199    psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
     200    psFree(filename);
     201    psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
     202    psFree(resolved);
     203    if (!fits || !pmReadoutReadSubtractionKernels(readoutCnv, fits)) {
     204        psError(PSPHOT_ERR_IO, false, "Unable to read previously produced kernel");
     205        psFitsClose(fits);
     206        return false;
     207    }
     208    psFitsClose(fits);
     209
     210    // read the convolved pixels (image, mask, variance) -- names are pre-defined
     211    if (!readImage(&readoutCnv->image,    options->convImages->data[index],    config) ||
     212        !readImage(&readoutCnv->mask,     options->convMasks->data[index],     config) ||
     213        !readImage(&readoutCnv->variance, options->convVariances->data[index], config)) {
     214        psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image.");
     215        return false;
     216    }
     217
     218    // XXX ??? not sure what is happening here -- consult Paul
     219    psRegion *region = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
     220    pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
     221
     222    pmSubtractionAnalysis(readoutCnv->analysis, NULL, kernels, region, readoutCnv->image->numCols, readoutCnv->image->numRows);
     223
     224    psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
     225
     226    // update the covariance matrix
     227    // XXX why is this needed if we have correctly read the saved data?
     228    bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
     229    psKernel *covar = psImageCovarianceCalculate(kernel, readoutCnv->covariance); // Covariance matrix
     230    psImageCovarianceSetThreads(oldThreads);
     231    psFree(readoutCnv->covariance);
     232    readoutCnv->covariance = covar;
     233    psFree(kernel);
     234    return true;
     235}
     236# endif
     237
     238bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname) {
     239
     240    pmHDU *hdu = pmHDUFromCell(readoutRef->parent);
     241    psString name = NULL;
     242    psStringAppend(&name, "%s_%03d.fits", rootname, index);
     243    pmStackVisualPlotTestImage(readoutOut->image, name);
     244    psFits *fits = psFitsOpen(name, "w");
     245    psFree(name);
     246    psFitsWriteImage(fits, hdu->header, readoutOut->image, 0, NULL);
     247    psFitsClose(fits);
     248    return true;
     249}
     250
     251bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname) {
     252
     253    pmHDU *hdu = pmHDUFromCell(readoutRef->parent);
     254    psString name = NULL;
     255    psStringAppend(&name, "%s_%03d.fits", rootname, index);
     256    pmStackVisualPlotTestImage(readoutFake->image, name);
     257    psFits *fits = psFitsOpen(name, "w");
     258    psFree(name);
     259    psBinaryOp(readoutFake->image, readoutConv->image, "-", readoutFake->image);
     260    psFitsWriteImage(fits, hdu->header, readoutFake->image, 0, NULL);
     261    psFitsClose(fits);
     262    return true;
     263}
     264
     265// perform the bulk of the PSF-matching
     266bool matchKernel(pmConfig *config, pmReadout *readoutOut, pmReadout *readoutSrc, psphotStackOptions *options, int index) {
     267
     268    bool mdok;
     269
     270    psAssert(options->psf, "Require target PSF");
     271    psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list");
     272
     273    psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
     274    psAssert(stackRecipe, "We've thrown an error on this before.");
     275
     276    // Look up appropriate values from the ppSub recipe
     277    psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
     278    psAssert(subRecipe, "recipe missing");
     279
     280    psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
     281    psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor
     282    psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
     283
     284    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     285    psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
     286    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
     287
     288    float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness
     289    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
     290
     291    int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order
     292    float regionSize = psMetadataLookupF32(NULL, subRecipe, "REGION.SIZE"); // Size of iso-kernel regs
     293    float spacing = psMetadataLookupF32(NULL, subRecipe, "STAMP.SPACING"); // Typical stamp spacing
     294
     295    int footprint = psMetadataLookupS32(NULL, subRecipe, "STAMP.FOOTPRINT"); // Stamp half-size
     296    int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size
     297
     298    float threshold = psMetadataLookupF32(NULL, subRecipe, "STAMP.THRESHOLD"); // Threshold for stmps
     299    int stride = psMetadataLookupS32(NULL, subRecipe, "STRIDE"); // Size of convolution patches
     300    int iter = psMetadataLookupS32(NULL, subRecipe, "ITER"); // Rejection iterations
     301    float rej = psMetadataLookupF32(NULL, subRecipe, "REJ"); // Rejection threshold
     302    float kernelError = psMetadataLookupF32(NULL, subRecipe, "KERNEL.ERR"); // Relative systematic error in kernel
     303    float normFrac = psMetadataLookupF32(NULL, subRecipe, "NORM.FRAC"); // Fraction of window for normalisn windw
     304    float sysError = psMetadataLookupF32(NULL, subRecipe, "SYS.ERR"); // Relative systematic error in images
     305    float skyErr = psMetadataLookupF32(NULL, subRecipe, "SKY.ERR"); // Additional error in sky
     306    float covarFrac = psMetadataLookupF32(NULL, subRecipe, "COVAR.FRAC"); // Fraction for covariance calculation
     307
     308    const char *typeStr = psMetadataLookupStr(NULL, subRecipe, "KERNEL.TYPE"); // Kernel type
     309    pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
     310    psVector *widths = psMetadataLookupPtr(NULL, subRecipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
     311    psVector *orders = psMetadataLookupPtr(NULL, subRecipe, "ISIS.ORDERS"); // ISIS Polynomial orders
     312    int inner = psMetadataLookupS32(NULL, subRecipe, "INNER"); // Inner radius
     313    int ringsOrder = psMetadataLookupS32(NULL, subRecipe, "RINGS.ORDER"); // RINGS polynomial order
     314    int binning = psMetadataLookupS32(NULL, subRecipe, "SPAM.BINNING"); // Binning for SPAM kernel
     315    float badFrac = psMetadataLookupF32(NULL, subRecipe, "BADFRAC"); // Maximum bad fraction
     316    float optThresh = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.TOL"); // Tolerance for search
     317    int optOrder = psMetadataLookupS32(&mdok, subRecipe, "OPTIMUM.ORDER"); // Order for search
     318    float poorFrac = psMetadataLookupF32(&mdok, subRecipe, "POOR.FRACTION"); // Fraction for "poor"
     319
     320    bool scale = psMetadataLookupBool(NULL, subRecipe, "SCALE");        // Scale kernel parameters?
     321    float scaleRef = psMetadataLookupF32(NULL, subRecipe, "SCALE.REF"); // Reference for scaling
     322    float scaleMin = psMetadataLookupF32(NULL, subRecipe, "SCALE.MIN"); // Minimum for scaling
     323    float scaleMax = psMetadataLookupF32(NULL, subRecipe, "SCALE.MAX"); // Maximum for scaling
     324    if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) {
     325        psError(PSPHOT_ERR_CONFIG, false,
     326                "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.",
     327                scaleRef, scaleMin, scaleMax);
     328        return false;
     329    }
     330
     331    // These values are specified specifically for stacking
     332    const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS");// Stamps filename
     333
     334    psVector *widthsCopy = NULL;
     335    psVector *optWidths = NULL;
     336    pmReadout *fake = NULL;
     337    psArray *stampSources = NULL;
     338
     339    bool optimum = false;
     340    optWidths = SetOptWidths(&optimum, subRecipe); // Vector with FWHMs for optimum search
     341
     342    // For the sake of stamps, remove nearby sources
     343    stampSources = stackSourcesFilter(options->sourceLists->data[index], footprint); // Filtered list of sources
     344
     345    fake = makeFakeReadout(config, readoutSrc, stampSources, options->psf, maskVal | maskBad, footprint + size);
     346    if (!fake) goto escape;
     347
     348    dumpImage(fake, readoutSrc, index, "fake");
     349    dumpImage(readoutSrc,  readoutSrc, index, "real");
     350
     351    if (threads) pmSubtractionThreadsInit();
     352
     353    // Do the image matching
     354    pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readoutSrc->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
     355    if (kernel) {
     356        if (!pmSubtractionMatchPrecalc(NULL, readoutOut, fake, readoutSrc, readoutSrc->analysis, stride, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac)) {
     357            psError(psErrorCodeLast(), false, "Unable to convolve images.");
     358            goto escape;
     359        }
     360    } else {
     361        // Scale the input parameters
     362        widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
     363        if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, options->inputSeeing->data.F32[index], options->targetSeeing, scaleRef, scaleMin, scaleMax)) {
     364            psError(psErrorCodeLast(), false, "Unable to scale kernel parameters");
     365            goto escape;
     366        }
     367
     368        if (!pmSubtractionMatch(NULL, readoutOut, fake, readoutSrc, footprint, stride, regionSize, spacing, threshold, stampSources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
     369            psError(psErrorCodeLast(), false, "Unable to match images.");
     370            goto escape;
     371        }
     372    }
     373
     374    // Reject image completely if the maximum deconvolution fraction exceeds the limit
     375    float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
     376    float deconv = psMetadataLookupF32(NULL, readoutOut->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
     377    if (deconv > deconvLimit) {
     378        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index);
     379        goto escape;
     380    }
     381
     382    dumpImage(readoutOut, readoutSrc, index, "conv");
     383    dumpImageDiff(readoutOut, fake, readoutSrc, index, "diff");
     384
     385    psFree(fake);
     386    psFree(optWidths);
     387    psFree(stampSources);
     388    psFree(widthsCopy);
     389    pmSubtractionThreadsFinalize();
     390    return true;
     391
     392escape:
     393    psFree(fake);
     394    psFree(optWidths);
     395    psFree(stampSources);
     396    psFree(widthsCopy);
     397    pmSubtractionThreadsFinalize();
     398    return false;
     399}
     400
     401// Extract the regions and solutions used in the image matching
     402// This stops them from being freed when we iterate back up the FPA
     403// Record the chi-square value
     404// XXX this function may not be needed for psphotStack
     405bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index) {
     406
     407    psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions
    287408    {
    288         pmHDU *hdu = pmHDUFromCell(readout->parent);
    289         psString name = NULL;
    290         psStringAppend(&name, "convolved_%03d.fits", index);
    291         pmStackVisualPlotTestImage(readout->image, name);
    292         psFits *fits = psFitsOpen(name, "w");
    293         psFree(name);
    294         psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
    295         psFitsClose(fits);
    296     }
    297 
    298 bool matchKernel() {
    299             // Normal operations here
    300             psAssert(options->psf, "Require target PSF");
    301             psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list");
    302 
    303             int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order
    304             float regionSize = psMetadataLookupF32(NULL, subRecipe, "REGION.SIZE"); // Size of iso-kernel regs
    305             float spacing = psMetadataLookupF32(NULL, subRecipe, "STAMP.SPACING"); // Typical stamp spacing
    306             int footprint = psMetadataLookupS32(NULL, subRecipe, "STAMP.FOOTPRINT"); // Stamp half-size
    307             float threshold = psMetadataLookupF32(NULL, subRecipe, "STAMP.THRESHOLD"); // Threshold for stmps
    308             int stride = psMetadataLookupS32(NULL, subRecipe, "STRIDE"); // Size of convolution patches
    309             int iter = psMetadataLookupS32(NULL, subRecipe, "ITER"); // Rejection iterations
    310             float rej = psMetadataLookupF32(NULL, subRecipe, "REJ"); // Rejection threshold
    311             float kernelError = psMetadataLookupF32(NULL, subRecipe, "KERNEL.ERR"); // Relative systematic error in kernel
    312             float normFrac = psMetadataLookupF32(NULL, subRecipe, "NORM.FRAC"); // Fraction of window for normalisn windw
    313             float sysError = psMetadataLookupF32(NULL, subRecipe, "SYS.ERR"); // Relative systematic error in images
    314             float skyErr = psMetadataLookupF32(NULL, subRecipe, "SKY.ERR"); // Additional error in sky
    315             float covarFrac = psMetadataLookupF32(NULL, subRecipe, "COVAR.FRAC"); // Fraction for covariance calculation
    316 
    317             const char *typeStr = psMetadataLookupStr(NULL, subRecipe, "KERNEL.TYPE"); // Kernel type
    318             pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
    319             psVector *widths = psMetadataLookupPtr(NULL, subRecipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
    320             psVector *orders = psMetadataLookupPtr(NULL, subRecipe, "ISIS.ORDERS"); // ISIS Polynomial orders
    321             int inner = psMetadataLookupS32(NULL, subRecipe, "INNER"); // Inner radius
    322             int ringsOrder = psMetadataLookupS32(NULL, subRecipe, "RINGS.ORDER"); // RINGS polynomial order
    323             int binning = psMetadataLookupS32(NULL, subRecipe, "SPAM.BINNING"); // Binning for SPAM kernel
    324             float badFrac = psMetadataLookupF32(NULL, subRecipe, "BADFRAC"); // Maximum bad fraction
    325             bool optimum = psMetadataLookupBool(&mdok, subRecipe, "OPTIMUM"); // Derive optimum parameters?
    326             float optMin = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.MIN"); // Minimum width for search
    327             float optMax = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.MAX"); // Maximum width for search
    328             float optStep = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.STEP"); // Step for search
    329             float optThresh = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.TOL"); // Tolerance for search
    330             int optOrder = psMetadataLookupS32(&mdok, subRecipe, "OPTIMUM.ORDER"); // Order for search
    331             float poorFrac = psMetadataLookupF32(&mdok, subRecipe, "POOR.FRACTION"); // Fraction for "poor"
    332 
    333             bool scale = psMetadataLookupBool(NULL, subRecipe, "SCALE");        // Scale kernel parameters?
    334             float scaleRef = psMetadataLookupF32(NULL, subRecipe, "SCALE.REF"); // Reference for scaling
    335             float scaleMin = psMetadataLookupF32(NULL, subRecipe, "SCALE.MIN"); // Minimum for scaling
    336             float scaleMax = psMetadataLookupF32(NULL, subRecipe, "SCALE.MAX"); // Maximum for scaling
    337             if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) {
    338                 psError(PPSTACK_ERR_CONFIG, false,
    339                         "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.",
    340                         scaleRef, scaleMin, scaleMax);
    341                 return false;
    342             }
    343 
    344 
    345             // These values are specified specifically for stacking
    346             const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename
    347 
    348             psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
    349             if (optimum) {
    350                 optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
    351             }
    352 
    353             pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
    354 
    355             psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background
    356             psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    357             if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
    358                 psError(PPSTACK_ERR_DATA, false, "Can't measure background for image.");
    359                 psFree(fake);
    360                 psFree(optWidths);
    361                 psFree(conv);
    362                 psFree(bg);
    363                 psFree(rng);
    364                 return false;
    365             }
    366             float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image
     409        psString regex = NULL;          // Regular expression
     410        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
     411        psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
     412        psFree(regex);
     413        psMetadataItem *item = NULL;// Item from iteration
     414        while ((item = psMetadataGetAndIncrement(iter))) {
     415            assert(item->type == PS_DATA_REGION);
     416            regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V);
     417        }
     418        psFree(iter);
     419    }
     420
     421    psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels
     422    {
     423        psString regex = NULL;          // Regular expression
     424        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
     425        psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
     426        psFree(regex);
     427        psMetadataItem *item = NULL;// Item from iteration
     428        while ((item = psMetadataGetAndIncrement(iter))) {
     429            assert(item->type == PS_DATA_UNKNOWN);
     430            pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
     431            kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
     432        }
     433        psFree(iter);
     434    }
     435    psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match");
     436
     437    // Record chi^2
     438    {
     439        double sum = 0.0;           // Sum of chi^2
     440        int num = 0;                // Number of measurements of chi^2
     441        psString regex = NULL;      // Regular expression
     442        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
     443        psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
     444        psFree(regex);
     445        psMetadataItem *item = NULL;// Item from iteration
     446        while ((item = psMetadataGetAndIncrement(iter))) {
     447            assert(item->type == PS_DATA_UNKNOWN);
     448            pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
     449            sum += kernels->mean;
     450            num++;
     451        }
     452        psFree(iter);
     453        options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
     454    }
     455
     456    return true;
     457}
     458
     459// Kernel normalisation for convolved readout
     460bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index) {
     461
     462    double sum = 0.0;           // Sum of chi^2
     463    int num = 0;                // Number of measurements of chi^2
     464    psString regex = NULL;      // Regular expression
     465    psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
     466    psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
     467    psFree(regex);
     468    psMetadataItem *item = NULL;// Item from iteration
     469    while ((item = psMetadataGetAndIncrement(iter))) {
     470        assert(item->type == PS_TYPE_F32);
     471        float norm = item->data.F32; // Normalisation
     472        sum += norm;
     473        num++;
     474    }
     475    psFree(iter);
     476    float conv = sum/num;       // Mean normalisation from convolution
     477    float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
     478    float renorm =  stars / conv; // Renormalisation to apply
     479    psLogMsg("psphotStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", index, renorm, conv, stars);
     480
     481    psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
     482    psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
     483    return true;
     484}
     485
     486// adjust scaling for readout (remove background, ..., determine weighting)
     487bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index) {
     488
     489    psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
     490    psAssert(stackRecipe, "We've thrown an error on this before.");
     491
     492    // Look up appropriate values from the ppSub recipe
     493    psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
     494    psAssert(subRecipe, "recipe missing");
     495
     496    psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
     497    psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
     498
     499    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     500    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
     501
     502    // Ensure the background value is zero
     503    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
     504    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     505
     506    // XXX why is this in config->arguments and not recipe?
     507    if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
     508        if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
     509            psAbort("Can't measure background for image.");
     510            // XXX we used to clear error: why is this acceptable? psErrorClear();
     511        }
     512
     513        float value = psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN);
     514        float stdev = psStatsGetValue(bg, PS_STAT_ROBUST_STDEV);
     515
     516        psLogMsg("psphotStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", value, stdev);
     517        psBinaryOp(readout->image, readout->image, "-", psScalarAlloc(value, PS_TYPE_F32));
     518    }
     519
     520    if (!stackRenormaliseReadout(config, readout)) {
     521        psFree(rng);
     522        psFree(bg);
     523        return false;
     524    }
     525
     526    // Measure the variance level for the weighting
     527    if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) {
     528        if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
     529            psError(PSPHOT_ERR_DATA, false, "Can't measure mean variance for image.");
    367530            psFree(rng);
    368531            psFree(bg);
    369 
    370             // For the sake of stamps, remove nearby sources
    371             psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index],
    372                                                        footprint); // Filtered list of sources
    373 
    374             bool oldThreads = pmReadoutFakeThreads(true); // Old threading state
    375             if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
    376                                           stampSources, SOURCE_MASK, NULL, NULL, options->psf,
    377                                           minFlux, footprint + size, false, true)) {
    378                 psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF.");
    379                 psFree(fake);
    380                 psFree(optWidths);
    381                 psFree(conv);
    382                 return false;
    383             }
    384             pmReadoutFakeThreads(oldThreads);
    385 
    386             fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK);
    387 
    388             // Add the background into the target image
    389             psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
    390             psBinaryOp(fake->image, fake->image, "+", bgImage);
    391             psFree(bgImage);
    392 
    393             dumpImage();
    394 
    395             if (threads > 0) {
    396                 pmSubtractionThreadsInit();
    397             }
    398 
    399             // Do the image matching
    400             pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readout->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
    401             if (kernel) {
    402                 if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis,
    403                                                stride, kernelError, covarFrac, maskVal, maskBad, maskPoor,
    404                                                poorFrac, badFrac)) {
    405                     psError(psErrorCodeLast(), false, "Unable to convolve images.");
    406                     psFree(fake);
    407                     psFree(optWidths);
    408                     psFree(stampSources);
    409                     psFree(conv);
    410                     if (threads > 0) {
    411                         pmSubtractionThreadsFinalize();
    412                     }
    413                     return false;
    414                 }
    415             } else {
    416                 // Scale the input parameters
    417                 psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
    418                 if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy,
    419                                                        options->inputSeeing->data.F32[index],
    420                                                        options->targetSeeing, scaleRef, scaleMin, scaleMax)) {
    421                     psError(psErrorCodeLast(), false, "Unable to scale kernel parameters");
    422                     psFree(fake);
    423                     psFree(optWidths);
    424                     psFree(stampSources);
    425                     psFree(conv);
    426                     psFree(widthsCopy);
    427                     if (threads > 0) {
    428                         pmSubtractionThreadsFinalize();
    429                     }
    430                     return false;
    431                 }
    432 
    433                 if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
    434                                         threshold, stampSources, stampsName, type, size, order, widthsCopy,
    435                                         orders, inner, ringsOrder, binning, penalty,
    436                                         optimum, optWidths, optOrder, optThresh, iter, rej, normFrac,
    437                                         sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor,
    438                                         poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
    439                     psError(psErrorCodeLast(), false, "Unable to match images.");
    440                     psFree(fake);
    441                     psFree(optWidths);
    442                     psFree(stampSources);
    443                     psFree(conv);
    444                     psFree(widthsCopy);
    445                     if (threads > 0) {
    446                         pmSubtractionThreadsFinalize();
    447                     }
    448                     return false;
    449                 }
    450                 psFree(widthsCopy);
    451             }
    452 
    453             dumpImage2();
    454 
    455             psFree(fake);
    456             psFree(optWidths);
    457             psFree(stampSources);
    458 
    459             if (threads > 0) {
    460                 pmSubtractionThreadsFinalize();
    461             }
    462 
    463             // Replace original images with convolved
    464             psFree(readout->image);
    465             psFree(readout->mask);
    466             psFree(readout->variance);
    467             psFree(readout->covariance);
    468             readout->image  = psMemIncrRefCounter(conv->image);
    469             readout->mask   = psMemIncrRefCounter(conv->mask);
    470             readout->variance = psMemIncrRefCounter(conv->variance);
    471             readout->covariance = psImageCovarianceTruncate(conv->covariance, COVAR_FRAC);
    472 
    473 }
    474 
    475 bool saveMatchData () {
    476        // Extract the regions and solutions used in the image matching
    477         // This stops them from being freed when we iterate back up the FPA
    478         psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions
    479         {
    480             psString regex = NULL;          // Regular expression
    481             psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
    482             psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
    483             psFree(regex);
    484             psMetadataItem *item = NULL;// Item from iteration
    485             while ((item = psMetadataGetAndIncrement(iter))) {
    486                 assert(item->type == PS_DATA_REGION);
    487                 regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V);
    488             }
    489             psFree(iter);
     532            return false;
    490533        }
    491         psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels
    492         {
    493             psString regex = NULL;          // Regular expression
    494             psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
    495             psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
    496             psFree(regex);
    497             psMetadataItem *item = NULL;// Item from iteration
    498             while ((item = psMetadataGetAndIncrement(iter))) {
    499                 assert(item->type == PS_DATA_UNKNOWN);
    500                 pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
    501                 kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
    502             }
    503             psFree(iter);
    504         }
    505         psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match");
    506 }
    507 
    508 bool saveChiSquare() {
    509         // Record chi^2
    510         {
    511             double sum = 0.0;           // Sum of chi^2
    512             int num = 0;                // Number of measurements of chi^2
    513             psString regex = NULL;      // Regular expression
    514             psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
    515             psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
    516             psFree(regex);
    517             psMetadataItem *item = NULL;// Item from iteration
    518             while ((item = psMetadataGetAndIncrement(iter))) {
    519                 assert(item->type == PS_DATA_UNKNOWN);
    520                 pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
    521                 sum += kernels->mean;
    522                 num++;
    523             }
    524             psFree(iter);
    525             options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
    526         }
    527 
    528 }
    529 
    530 bool renormKernel() {
    531         // Kernel normalisation
    532         {
    533             double sum = 0.0;           // Sum of chi^2
    534             int num = 0;                // Number of measurements of chi^2
    535             psString regex = NULL;      // Regular expression
    536             psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
    537             psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
    538             psFree(regex);
    539             psMetadataItem *item = NULL;// Item from iteration
    540             while ((item = psMetadataGetAndIncrement(iter))) {
    541                 assert(item->type == PS_TYPE_F32);
    542                 float norm = item->data.F32; // Normalisation
    543                 sum += norm;
    544                 num++;
    545             }
    546             psFree(iter);
    547             float conv = sum/num;       // Mean normalisation from convolution
    548             float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
    549             float renorm =  stars / conv; // Renormalisation to apply
    550             psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n",
    551                      index, renorm, conv, stars);
    552             psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
    553             psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
    554         }
    555 
    556 }
     534        options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance));
     535    } else {
     536        options->weightings->data.F32[index] = 1.0;
     537    }
     538    psLogMsg("psphotStack", PS_LOG_INFO, "Weighting for image %d is %f\n", index, options->weightings->data.F32[index]);
     539
     540    psFree(rng);
     541    psFree(bg);
     542    return true;
     543}
     544
     545# define NOISE_FRACTION 0.01             // Set minimum flux to this fraction of noise
     546# define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
     547
     548// generate a fake readout against which to PSF match
     549pmReadout *makeFakeReadout(pmConfig *config, pmReadout *readoutSrc, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize) {
     550
     551    pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
     552
     553    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background
     554    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     555    if (!psImageBackground(bg, NULL, readoutSrc->image, readoutSrc->mask, maskVal, rng)) {
     556        psError(PSPHOT_ERR_DATA, false, "Can't measure background for image.");
     557        psFree(fake);
     558        psFree(bg);
     559        psFree(rng);
     560        return NULL;
     561    }
     562    float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image
     563    psFree(rng);
     564    psFree(bg);
     565
     566    bool oldThreads = pmReadoutFakeThreads(true); // Old threading state
     567    if (!pmReadoutFakeFromSources(fake, readoutSrc->image->numCols, readoutSrc->image->numRows, sources, SOURCE_MASK, NULL, NULL, psf, minFlux, fullSize, false, true)) {
     568        psError(PSPHOT_ERR_DATA, false, "Unable to generate fake image with target PSF.");
     569        psFree(fake);
     570        return NULL;
     571    }
     572    pmReadoutFakeThreads(oldThreads);
     573
     574    fake->mask = psImageCopy(NULL, readoutSrc->mask, PS_TYPE_IMAGE_MASK);
     575
     576    // Add the background into the target image
     577    psImage *bgImage = stackBackgroundModel(readoutSrc, config); // Image of background
     578    psBinaryOp(fake->image, fake->image, "+", bgImage);
     579    psFree(bgImage);
     580
     581    return fake;
     582}
     583
     584// set the widths
     585psVector *SetOptWidths (bool *optimum, psMetadata *recipe) {
     586
     587    bool status;
     588
     589    *optimum = psMetadataLookupBool(&status, recipe, "OPTIMUM"); // Derive optimum parameters?
     590    psAssert (status, "missing recipe value %s", "OPTIMUM");
     591
     592    psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
     593
     594    if (*optimum) {
     595        float optMin = psMetadataLookupF32(&status,  recipe, "OPTIMUM.MIN"); // Minimum width for search
     596        psAssert (status, "missing recipe value %s", "OPTIMUM.MIN");
     597       
     598        float optMax = psMetadataLookupF32(&status,  recipe, "OPTIMUM.MAX"); // Maximum width for search
     599        psAssert (status, "missing recipe value %s", "OPTIMUM.MAX");
     600       
     601        float optStep = psMetadataLookupF32(&status, recipe, "OPTIMUM.STEP"); // Step for search
     602        psAssert (status, "missing recipe value %s", "OPTIMUM.STEP");
     603
     604        optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
     605    }
     606
     607    return optWidths;
     608}

  • trunk/psphot/src/psphotStackParseCamera.c

    r27657 r28013  
    2525        psMetadata *input = item->data.md; // The input metadata of interest
    2626
    27         // look for 'IMAGE', 'MASK', 'VARIANCE' in folder (only 'IMAGE' is required)
    28 
    29         psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image
    30         if (!image || strlen(image) == 0) {
    31             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name);
     27        pmFPAfile *rawInputFile = NULL;
     28        pmFPAfile *cnvInputFile = NULL;
     29
     30        // RAW (unconvolved) input data (RAW:IMAGE, RAW:MASK, RAW:VARIANCE, RAW:PSF)
     31        psString rawImage = psMetadataLookupStr(&status, input, "RAW:IMAGE"); // Name of image
     32        if (rawImage && strlen(rawImage) > 0) {
     33            rawInputFile = defineFile(config, NULL, "PSPHOT.STACK.INPUT.RAW", rawImage, PM_FPA_FILE_IMAGE); // File for image
     34            if (!rawInputFile) {
     35                psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, rawImage);
     36                return false;
     37            }
     38            psString mask = psMetadataLookupStr(&status, input, "RAW:MASK"); // Name of mask
     39            if (mask && strlen(mask) > 0) {
     40                if (!defineFile(config, rawInputFile, "PSPHOT.STACK.MASK.RAW", mask, PM_FPA_FILE_MASK)) {
     41                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
     42                    return false;
     43                }
     44            }
     45            psString variance = psMetadataLookupStr(&status, input, "RAW:VARIANCE"); // Name of variance map
     46            if (variance && strlen(variance) > 0) {
     47                if (!defineFile(config, rawInputFile, "PSPHOT.STACK.VARIANCE.RAW", variance, PM_FPA_FILE_VARIANCE)) {
     48                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance);
     49                    return false;
     50                }
     51            }
     52            psString psf = psMetadataLookupStr(&status, input, "RAW:PSF"); // Name of mask
     53            if (psf && strlen(psf) > 0) {
     54                if (!defineFile(config, rawInputFile, "PSPHOT.STACK.PSF.RAW", psf, PM_FPA_FILE_PSF)) {
     55                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf);
     56                    return false;
     57                }
     58            }
     59        }
     60
     61        // CNV (convolved) input data (CNV:IMAGE, CNV:MASK, CNV:VARIANCE, CNV:PSF)
     62        psString cnvImage = psMetadataLookupStr(&status, input, "CNV:IMAGE"); // Name of image
     63        if (cnvImage && strlen(cnvImage) > 0) {
     64            cnvInputFile = defineFile(config, NULL, "PSPHOT.STACK.INPUT.CNV", cnvImage, PM_FPA_FILE_IMAGE); // File for image
     65            if (!cnvInputFile) {
     66                psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, cnvImage);
     67                return false;
     68            }
     69            psString mask = psMetadataLookupStr(&status, input, "CNV:MASK"); // Name of mask
     70            if (mask && strlen(mask) > 0) {
     71                if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.MASK.CNV", mask, PM_FPA_FILE_MASK)) {
     72                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
     73                    return false;
     74                }
     75            }
     76            psString variance = psMetadataLookupStr(&status, input, "CNV:VARIANCE"); // Name of variance map
     77            if (variance && strlen(variance) > 0) {
     78                if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.VARIANCE.CNV", variance, PM_FPA_FILE_VARIANCE)) {
     79                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance);
     80                    return false;
     81                }
     82            }
     83            psString psf = psMetadataLookupStr(&status, input, "CNV:PSF"); // Name of mask
     84            if (psf && strlen(psf) > 0) {
     85                if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.PSF.CNV", psf, PM_FPA_FILE_PSF)) {
     86                    psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf);
     87                    return false;
     88                }
     89            }
     90        }
     91
     92        if (!rawInputFile && !cnvInputFile) {
     93            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s (%d) lacks IMAGE of type STR", item->name, i);
    3294            return false;
    3395        }
    34         pmFPAfile *imageFile = defineFile(config, NULL, "PSPHOT.INPUT", image, PM_FPA_FILE_IMAGE); // File for image
    35         if (!imageFile) {
    36             psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image);
    37             return false;
    38         }
    39 
    40         psString mask = psMetadataLookupStr(&status, input, "MASK"); // Name of mask
    41         if (mask && strlen(mask) > 0) {
    42             if (!defineFile(config, imageFile, "PSPHOT.MASK", mask, PM_FPA_FILE_MASK)) {
    43                 psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
    44                 return false;
    45             }
    46         }
    47 
    48         psString variance = psMetadataLookupStr(&status, input, "VARIANCE"); // Name of variance map
    49         if (variance && strlen(variance) > 0) {
    50             if (!defineFile(config, imageFile, "PSPHOT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) {
    51                 psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance);
    52                 return false;
    53             }
    54         }
    55         // the output sources are carried on the input->fpa structures
    56         pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, imageFile, "PSPHOT.STACK.OUTPUT");
    57         if (!outsources) {
    58             psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT");
    59             return false;
    60         }
    61         outsources->save = true;
    62         outsources->fileID = i;         // this is used to generate output names
     96
     97        psString sources = psMetadataLookupStr(&status, input, "SOURCES"); // Name of mask
     98        // pmFPAfile *srcInputFile = rawInputFile ? rawInputFile : cnvInputFile;
     99        if (sources && strlen(sources) > 0) {
     100            if (!defineFile(config, NULL, "PSPHOT.STACK.SOURCES", sources, PM_FPA_FILE_CMF)) {
     101                psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)", i, sources);
     102                return false;
     103            }
     104        }
     105
     106        // generate an pmFPAimage for the output convolved image
     107        // XXX output of these files should be optional
     108        {
     109            pmFPAfile *outputImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.STACK.OUTPUT.IMAGE");
     110            if (!outputImage) {
     111                psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.STACK.OUTPUT.IMAGE");
     112                return false;
     113            }
     114            outputImage->save = true;
     115            outputImage->fileID = i;            // this is used to generate output names
     116
     117            pmFPAfile *outputMask = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.MASK");
     118            if (!outputMask) {
     119                psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.STACK.OUTPUT.MASK"));
     120                return NULL;
     121            }
     122            if (outputMask->type != PM_FPA_FILE_MASK) {
     123                psError(PS_ERR_IO, true, "PSPHOT.STACK.OUTPUT.MASK is not of type MASK");
     124                return NULL;
     125            }
     126            outputMask->save = true;
     127            outputMask->fileID = i;             // this is used to generate output names
     128
     129            pmFPAfile *outputVariance = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.VARIANCE");
     130            if (!outputVariance) {
     131                psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.STACK.OUTPUT.VARIANCE"));
     132                return NULL;
     133            }
     134            if (outputVariance->type != PM_FPA_FILE_VARIANCE) {
     135                psError(PS_ERR_IO, true, "PSPHOT.STACK.OUTPUT.VARIANCE is not of type VARIANCE");
     136                return NULL;
     137            }
     138            outputVariance->save = true;
     139            outputVariance->fileID = i;         // this is used to generate output names
     140
     141            // the output sources are carried on the outputImage->fpa structures
     142            pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, outputImage, "PSPHOT.STACK.OUTPUT");
     143            if (!outsources) {
     144                psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT");
     145                return false;
     146            }
     147            outsources->save = true;
     148            outsources->fileID = i;             // this is used to generate output names
     149        }
    63150    }
    64151    psMetadataRemoveKey(config->arguments, "FILENAMES");
    … …  
    71158
    72159    // generate an pmFPAimage for the chisqImage
    73     pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE");
    74     if (!chisqImage) {
    75         psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE");
    76         return false;
    77     }
    78     chisqImage->save = true;
    79 
    80     pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK");
    81     if (!chisqMask) {
    82         psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK"));
    83         return NULL;
    84     }
    85     if (chisqMask->type != PM_FPA_FILE_MASK) {
    86         psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK");
    87         return NULL;
    88     }
    89     chisqMask->save = true;
    90 
    91     pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE");
    92     if (!chisqVariance) {
    93         psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE"));
    94         return NULL;
    95     }
    96     if (chisqVariance->type != PM_FPA_FILE_VARIANCE) {
    97         psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE");
    98         return NULL;
    99     }
    100     chisqVariance->save = true;
    101 
    102 # if (0)   
    103     // define the additional input/output files associated with psphot
    104     // XXX figure out which files are needed by psphotStack
    105     if (false && !psphotDefineFiles (config, input)) {
    106         psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files");
    107         return false;
    108     }
    109 # endif
     160    // XXX output of these files should be optional
     161    {
     162        pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE");
     163        if (!chisqImage) {
     164            psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE");
     165            return false;
     166        }
     167        chisqImage->save = true;
     168
     169        pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK");
     170        if (!chisqMask) {
     171            psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK"));
     172            return NULL;
     173        }
     174        if (chisqMask->type != PM_FPA_FILE_MASK) {
     175            psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK");
     176            return NULL;
     177        }
     178        chisqMask->save = true;
     179
     180        pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE");
     181        if (!chisqVariance) {
     182            psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE"));
     183            return NULL;
     184        }
     185        if (chisqVariance->type != PM_FPA_FILE_VARIANCE) {
     186            psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE");
     187            return NULL;
     188        }
     189        chisqVariance->save = true;
     190    }
    110191
    111192    psTrace("psphot", 1, "Done with psphotStackParseCamera...\n");
    … …  
    145226}
    146227
     228
     229/***
     230 *
     231 *  psphotStack :
     232
     233 *    * inputs:
     234 *      * unconvolved images
     235 *      * raw convolved images
     236 *      * psfs (unconvolved or convolved?)
     237 *      * sources
     238 
     239 * optionally convolve the unconvolved or the raw inputs
     240 * optionally perform no convolutions
     241 * optionally save the psf-matched images
     242
     243 */
     244
     245   
     246# if (0)   
     247    // define the additional input/output files associated with psphot
     248    // XXX figure out which files are needed by psphotStack
     249    if (false && !psphotDefineFiles (config, input)) {
     250        psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files");
     251        return false;
     252    }
     253# endif
     254

  • trunk/psphot/src/psphotStackReadout.c

    r27657 r28013  
    11# include "psphotInternal.h"
     2
     3# define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
     4# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"
    25
    36bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
    … …  
    1720    PS_ASSERT_PTR_NON_NULL (breakPt, false);
    1821
     22    // we have 3 relevant files: RAW, CNV, OUT
     23
    1924    // set the photcode for each image
    20     if (!psphotAddPhotcode (config, view)) {
     25    if (!psphotAddPhotcode (config, view, STACK_OUT)) {
    2126        psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
    2227        return false;
    … …  
    2429
    2530    // Generate the mask and weight images
    26     if (!psphotSetMaskAndVariance (config, view)) {
    27         return psphotReadoutCleanup (config, view);
     31    // XXX this should be done before we perform the convolutions
     32    if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) {
     33        return psphotReadoutCleanup (config, view, STACK_OUT);
    2834    }
    2935    if (!strcasecmp (breakPt, "NOTHING")) {
    30         return psphotReadoutCleanup (config, view);
     36        return psphotReadoutCleanup (config, view, STACK_OUT);
    3137    }
    3238
    3339    // generate a background model (median, smoothed image)
    3440    // XXX I think this is not defined correctly for an array of images.
    35     if (!psphotModelBackground (config, view)) {
    36         return psphotReadoutCleanup (config, view);
     41    // XXX probably need to subtract the model (same model?) for both RAW and OUT
     42    if (!psphotModelBackground (config, view, STACK_RAW)) {
     43        return psphotReadoutCleanup (config, view, STACK_OUT);
    3744    }
    38     if (!psphotSubtractBackground (config, view)) {
    39         return psphotReadoutCleanup (config, view);
     45    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
     46        return psphotReadoutCleanup (config, view, STACK_OUT);
    4047    }
    4148    if (!strcasecmp (breakPt, "BACKMDL")) {
    42         return psphotReadoutCleanup (config, view);
     49        return psphotReadoutCleanup (config, view, STACK_OUT);
    4350    }
    4451
    … …  
    4754    if (!psphotLoadPSF (config, view)) {
    4855        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
    49         return psphotReadoutCleanup (config, view);
     56        return psphotReadoutCleanup (config, view, STACK_OUT);
    5057    }
    5158
    52     if (!psphotStackChisqImage(config, view)) {
     59    if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_OUT)) {
    5360        psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
    54         return psphotReadoutCleanup (config, view);
     61        return psphotReadoutCleanup (config, view, STACK_OUT);
    5562    }
    5663    if (!strcasecmp (breakPt, "CHISQ")) {
    57         return psphotReadoutCleanup (config, view);
     64        return psphotReadoutCleanup (config, view, STACK_OUT);
    5865    }
    5966
    6067    // find the detections (by peak and/or footprint) in the image.
    6168    // This finds the detections on Chisq image as well as the individuals
    62     if (!psphotFindDetections (config, view, true)) { // pass 1
     69    if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
    6370        // this only happens if we had an error in psphotFindDetections
    6471        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    65         return psphotReadoutCleanup (config, view);
     72        return psphotReadoutCleanup (config, view, STACK_OUT);
     73    }
     74
     75    // copy the detections from RAW to OUT
     76    if (!psphotCopySources (config, view, STACK_OUT, STACK_RAW)) {
     77        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
     78        return psphotReadoutCleanup (config, view, STACK_OUT);
    6679    }
    6780
    6881    // construct sources and measure basic stats (saved on detections->newSources)
    6982    // only run this on detections from the input images, not chisq image
    70     if (!psphotSourceStats (config, view, true)) { // pass 1
     83    if (!psphotSourceStats (config, view, STACK_OUT, true)) { // pass 1
    7184        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    72         return psphotReadoutCleanup (config, view);
     85        return psphotReadoutCleanup (config, view, STACK_OUT);
    7386    }
    7487
    75     // *** generate the objects (which unify the sources from the different images)
    76     psArray *objects = psphotMatchSources (config, view);
     88    // generate the objects (object unify the sources from the different images)
     89    psArray *objects = psphotMatchSources (config, view, STACK_OUT);
    7790
    7891    // construct sources for the newly-generated sources (from other images)
    79     if (!psphotSourceStats (config, view, false)) { // pass 1
     92    if (!psphotSourceStats (config, view, STACK_OUT, false)) { // pass 1
    8093        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    81         return psphotReadoutCleanup (config, view);
     94        return psphotReadoutCleanup (config, view, STACK_OUT);
    8295    }
    8396
    … …  
    8598    // if (!psphotDeblendSatstars (config, view)) {
    8699    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
    87     //     return psphotReadoutCleanup (config, view);
     100    //     return psphotReadoutCleanup (config, view, STACK_OUT);
    88101    // }
    89102
    … …  
    91104    // if (!psphotBasicDeblend (config, view)) {
    92105    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
    93     //     return psphotReadoutCleanup (config, view);
     106    //     return psphotReadoutCleanup (config, view, STACK_OUT);
    94107    // }
    95108
    96109    // classify sources based on moments, brightness
    97110    // only run this on detections from the input images, not chisq image
    98     if (!psphotRoughClass (config, view)) {
     111    if (!psphotRoughClass (config, view, STACK_OUT)) {
    99112        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
    100         return psphotReadoutCleanup (config, view);
     113        return psphotReadoutCleanup (config, view, STACK_OUT);
    101114    }
    102115    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
    103116    // only run this on detections from the input images, not chisq image
    104     if (!psphotImageQuality (config, view)) { // pass 1
     117    if (!psphotImageQuality (config, view, STACK_OUT)) { // pass 1
    105118        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
    106         return psphotReadoutCleanup(config, view);
     119        return psphotReadoutCleanup (config, view, STACK_OUT);
    107120    }
    108121    if (!strcasecmp (breakPt, "MOMENTS")) {
    109         return psphotReadoutCleanup (config, view);
     122        return psphotReadoutCleanup (config, view, STACK_OUT);
    110123    }
    111124
    112125    // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
    113126    // this step is skipped
    114     if (!psphotChoosePSF (config, view)) { // pass 1
     127    if (!psphotChoosePSF (config, view, STACK_OUT)) { // pass 1
    115128        psLogMsg ("psphot", 3, "failure to construct a psf model");
    116         return psphotReadoutCleanup (config, view);
     129        return psphotReadoutCleanup (config, view, STACK_OUT);
    117130    }
    118131    if (!strcasecmp (breakPt, "PSFMODEL")) {
    119         return psphotReadoutCleanup (config, view);
     132        return psphotReadoutCleanup (config, view, STACK_OUT);
    120133    }
    121134
    122     // include externally-supplied sources
    123     // XXX fix this in the new multi-input context
    124     // psphotLoadExtSources (config, view); // pass 1
    125 
    126135    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    127     psphotGuessModels (config, view);
     136    psphotGuessModels (config, view, STACK_OUT);
    128137
    129138    // merge the newly selected sources into the existing list
    130139    // NOTE: merge OLD and NEW
    131     psphotMergeSources (config, view);
     140    psphotMergeSources (config, view, STACK_OUT);
    132141
    133142    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    134143    psphotFitSourcesLinearStack (config, objects, FALSE);
    135     psFree (objects);
    136144
    137145    // identify CRs and extended sources
    138     psphotSourceSize (config, view, TRUE);
     146    psphotSourceSize (config, view, STACK_OUT, TRUE);
    139147
    140148    // measure aperture photometry corrections
    141     if (!psphotApResid (config, view)) {
     149    if (!psphotApResid (config, view, STACK_OUT)) {
     150        psFree (objects);
    142151        psLogMsg ("psphot", 3, "failed on psphotApResid");
    143         return psphotReadoutCleanup (config, view);
     152        return psphotReadoutCleanup (config, view, STACK_OUT);
    144153    }
    145154
     155    psphotStackObjectsUnifyPosition (objects);
     156    psphotRadialAperturesByObject (config, objects, view, STACK_OUT);
     157
     158    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_OUT); // pass 1 (detections->allSources)
     159    psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources)
     160
    146161    // calculate source magnitudes
    147     psphotMagnitudes(config, view);
     162    psphotMagnitudes(config, view, STACK_OUT);
    148163
    149     if (!psphotEfficiency(config, view)) {
     164    if (!psphotEfficiency(config, view, STACK_OUT)) {
    150165        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    151166        psErrorClear();
    … …  
    156171
    157172    // replace background in residual image
    158     psphotSkyReplace (config, view);
     173    psphotSkyReplace (config, view, STACK_RAW);
    159174
    160175    // drop the references to the image pixels held by each source
    161     psphotSourceFreePixels (config, view);
     176    psphotSourceFreePixels (config, view, STACK_OUT);
    162177
    163178    // remove chisq image from config->file:PSPHOT.INPUT (why?)
    164     psphotStackRemoveChisqFromInputs(config);
     179    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
     180
     181    psFree (objects);
    165182
    166183    // create the exported-metadata and free local data
    167     return psphotReadoutCleanup (config, view);
     184    return psphotReadoutCleanup (config, view, STACK_OUT);
    168185}
    169186

  • trunk/psphot/src/psphotSubtractBackground.c

    r27657 r28013  
    44// generate the median in NxN boxes, clipping heavily
    55// linear interpolation to generate full-scale model
    6 bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
     6bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe)
    77{
    88    bool status = true;
    … …  
    1313
    1414    // find the currently selected readout
    15     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     15    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    1616
    1717    pmFPA *inFPA = file->fpa;
    … …  
    124124}
    125125
    126 bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view)
     126bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filerule)
    127127{
    128128    bool status = false;
    … …  
    137137    // loop over the available readouts
    138138    for (int i = 0; i < num; i++) {
    139         if (!psphotSubtractBackgroundReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     139        if (!psphotSubtractBackgroundReadout (config, view, filerule, i, recipe)) {
    140140            psError (PSPHOT_ERR_CONFIG, false, "failed to subtract background for PSPHOT.INPUT entry %d", i);
    141141            return false;
Note: See TracChangeset for help on using the changeset viewer.