IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30101


Ignore:
Timestamp:
Dec 17, 2010, 9:57:34 AM (15 years ago)
Author:
eugene
Message:

various modifications to psphotStack: add needed metadata to header; enable multiple PSF-matched images; work to speed up radial aperture photometry; radialAper is now an array (one element per PSF size); better feedback on extended fit skips; clarify processing: main analysis is on the input images (raw or convolved) while radial apertures is on the psf-matched images

Location:
branches/eam_branches/ipp-20101205/psphot/src
Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20101205/psphot/src/Makefile.am

    r29936 r30101  
    9898        psphotStackMatchPSFsUtils.c   \
    9999        psphotStackMatchPSFsPrepare.c \
     100        psphotStackMatchPSFsNext.c    \
    100101        psphotStackOptions.c          \
    101102        psphotStackObjects.c          \

  • branches/eam_branches/ipp-20101205/psphot/src/psphot.h

    r29936 r30101  
    312312
    313313int psphotFileruleCount(const pmConfig *config, const char *filerule);
     314bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num);
    314315
    315316bool psphotAddKnownSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *inSources);
    … …  
    403404bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
    404405bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
     406bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int nextSize);
     407bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, float currentFWHM, float targetFWHM);
    405408
    406409// psphotStackMatchPSFsUtils
    … …  
    426429bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule);
    427430bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    428 bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax);
     431bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry);
    429432
    430433bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
    431 bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
     434bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule, int nMatchedPSF);
    432435
    433436bool psphotStackObjectsUnifyPosition (psArray *objects);

  • branches/eam_branches/ipp-20101205/psphot/src/psphotExtendedSourceAnalysisByObject.c

    r29936 r30101  
    5353        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    5454        psAssert (readout, "missing readout?");
     55
     56        psLogMsg("psphot", PS_LOG_INFO, "petrosians for image %d", i);
     57        psphotVisualShowImage(readout);
    5558
    5659        readouts->data[i] = psMemIncrRefCounter(readout);

  • branches/eam_branches/ipp-20101205/psphot/src/psphotExtendedSourceFits.c

    r29936 r30101  
    1818    int num = psphotFileruleCount(config, filerule);
    1919
     20    // skip the chisq image (optionally?)
     21    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     22    if (!status) chisqNum = -1;
     23
    2024    // loop over the available readouts
    2125    for (int i = 0; i < num; i++) {
     26        if (i == chisqNum) continue; // skip chisq image
    2227        if (!psphotExtendedSourceFitsReadout (config, view, filerule, i, recipe)) {
    2328            psError (PSPHOT_ERR_CONFIG, false, "failed on to fit extended sources for %s entry %d", filerule, i);
    … …  
    3742    int Nplain = 0;
    3843    int NplainPass = 0;
     44    int Nfaint = 0;
    3945
    4046    psTimerStart ("psphot.extended");
    … …  
    4652    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    4753    psAssert (readout, "missing readout?");
     54
     55    psLogMsg("psphot", PS_LOG_INFO, "extended source fits for image %d", index);
     56    psphotVisualShowImage(readout);
    4857
    4958    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    … …  
    167176            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nplain
    168177            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass
     178            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfain
    169179
    170180            if (false && !psThreadJobAddPending(job)) {
    … …  
    189199                scalar = job->args->data[11];
    190200                NplainPass += scalar->data.S32;
     201                scalar = job->args->data[12];
     202                Nfaint += scalar->data.S32;
    191203                psFree(job);
    192204            }
    … …  
    217229                scalar = job->args->data[11];
    218230                NplainPass += scalar->data.S32;
     231                scalar = job->args->data[12];
     232                Nfaint += scalar->data.S32;
    219233            }
    220234            psFree(job);
    … …  
    227241    psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
    228242    psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
     243    psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit\n", Nfaint);
    229244    return true;
    230245}
    … …  
    238253    int NconvolvePass = 0;
    239254    int Nplain = 0;
     255    int Nfaint = 0;
    240256    int NplainPass = 0;
    241257    bool savePics = false;
    … …  
    350366          // limit selection to some SN limit
    351367          assert (source->peak); // how can a source not have a peak?
    352           if (source->peak->SN < SNlim) continue;
     368          if (source->peak->SN < SNlim) {
     369              Nfaint ++;
     370              continue;
     371          }
    353372
    354373          // check on the model type
    … …  
    497516    scalar->data.S32 = NplainPass;
    498517
     518    scalar = job->args->data[12];
     519    scalar->data.S32 = Nfaint;
     520
    499521    return true;
    500522}

  • branches/eam_branches/ipp-20101205/psphot/src/psphotRadialApertures.c

    r29936 r30101  
    102102        if (source->peak->x > AnalysisRegion.x1) continue;
    103103        if (source->peak->y > AnalysisRegion.y1) continue;
     104
     105        // allocate pmSourceExtendedParameters, if not already defined
     106        if (!source->radialAper) {
     107            source->radialAper = psArrayAlloc(1);
     108        }
    104109
    105110        // replace object in image
    … …  
    116121        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
    117122
    118         if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) {
     123        if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, 0)) {
    119124            psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    120125        } else {
    … …  
    130135}
    131136
    132 bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax) {
    133 
    134     // allocate pmSourceExtendedParameters, if not already defined
    135     if (!source->radial) {
    136         source->radial = pmSourceRadialAperturesAlloc ();
     137bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
     138
     139    psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?");
     140
     141    pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc ();
     142    source->radialAper->data[entry] = radialAper;
     143
     144    // storage for the derived pixel values
     145    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
     146    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
     147    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
     148
     149    // outer-most radius for initial truncation
     150    float Rmax  = aperRadii->data.F32[aperRadii->n - 1];
     151    float Rmax2 = PS_SQR(Rmax);
     152
     153    // store the R^2 values for the apertures
     154    psVector *aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
     155    for (int i = 0; i < aperRadii->n; i++) {
     156        aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
     157    }
     158
     159    // center of the apertures
     160    float xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
     161    float yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
     162
     163    // one pass through the pixels to select the valid pixels and calculate R^2
     164    for (int iy = 0; iy < source->pixels->numRows; iy++) {
     165
     166        float yDiff = iy - yCM;
     167        if (fabs(yDiff) > Rmax) continue;
     168
     169        float *vPix = source->pixels->data.F32[iy];
     170        float *vWgt = source->variance->data.F32[iy];
     171        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy];
     172
     173        for (int ix = 0; ix < source->pixels->numCols; ix++, vPix++, vWgt++) {
     174
     175            if (vMsk) {
     176                if (*vMsk & maskVal) {
     177                    vMsk++;
     178                    continue;
     179                }
     180                vMsk++;
     181            }
     182            if (isnan(*vPix)) continue;
     183
     184            float xDiff = ix - xCM;
     185            if (fabs(xDiff) > Rmax) continue;
     186
     187            // radius is just a function of (xDiff, yDiff)
     188            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     189            if (r2 > Rmax2) continue;
     190
     191            psVectorAppend(pixRadius2, r2);
     192            psVectorAppend(pixFlux, *vPix);
     193            psVectorAppend(pixVar, *vWgt);
     194        }
     195    }
     196
     197    psVector *flux    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
     198    psVector *fluxErr = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
     199    psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
     200
     201    psVectorInit (flux,    0.0);
     202    psVectorInit (fluxErr, 0.0);
     203    psVectorInit (fill,    0.0);
     204
     205    float *rPix2 = pixRadius2->data.F32;
     206    for (int i = 0; i < pixRadius2->n; i++, rPix2++) {
     207
     208        float *aRad2 = aperRadii2->data.F32;
     209        for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) {
     210            flux->data.F32[j]    += pixFlux->data.F32[i];
     211            fluxErr->data.F32[j] += pixVar->data.F32[i];
     212            fill->data.F32[j]    += 1.0;
     213        }
     214    }
     215
     216    for (int i = 0; i < flux->n; i++) {
     217        // calculate the total flux for bin 'nOut'
     218        float Area = M_PI*aperRadii2->data.F32[i];
     219        fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
     220        fill->data.F32[i] /= Area;
     221        psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n",
     222                 i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], fill->data.F32[i], Area);
    137223    }
    138224   
    139     psVector *radius  = psVectorAllocEmpty(100, PS_TYPE_F32);
     225    radialAper->flux = flux;
     226    radialAper->fluxErr = fluxErr;
     227    radialAper->fill = fill;
     228
     229    psFree (aperRadii2);
     230    psFree (pixRadius2);
     231    psFree (pixFlux);
     232    psFree (pixVar);
     233
     234    return true;
     235}
     236
     237static int nCalls = 0;
     238static int nPass = 0;
     239static int nPix = 0;
     240
     241bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) {
     242
     243    psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?");
     244
     245    pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc ();
     246    source->radialAper->data[entry] = radialAper;
     247
     248    psVector *pixRadius  = psVectorAllocEmpty(100, PS_TYPE_F32);
    140249    psVector *pixFlux = psVectorAllocEmpty(100, PS_TYPE_F32);
    141250    psVector *pixVar  = psVectorAllocEmpty(100, PS_TYPE_F32);
    … …  
    144253        for (int ix = 0; ix < source->pixels->numCols; ix++) {
    145254
    146             // 0.5 PIX: get radius as a function of pixel coord
     255            // 0.5 PIX: get pixRadius as a function of pixel coord
    147256            float x = ix + 0.5 - source->peak->xf + source->pixels->col0;
    148257            float y = iy + 0.5 - source->peak->yf + source->pixels->row0;
    … …  
    150259            float r = hypot(x, y);
    151260
    152             psVectorAppend(radius, r);
     261            psVectorAppend(pixRadius, r);
    153262            psVectorAppend(pixFlux, source->pixels->data.F32[iy][ix]);
    154263            psVectorAppend(pixVar, source->variance->data.F32[iy][ix]);
    155         }
    156     }
    157     psphotRadialAperturesSortFlux(radius, pixFlux, pixVar);
     264            nPix ++;
     265            // if (nPix % 10000 == 0) {fprintf (stderr, "?");}
     266        }
     267    }
     268    psphotRadialAperturesSortFlux(pixRadius, pixFlux, pixVar);
    158269
    159270    psVector *flux    = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin
    … …  
    174285
    175286    // XXX assume (or enforce) that the bins are contiguous and non-overlapping (Rmax[i] = Rmin[i+1])
    176     for (int i = 0; !done && (i < radius->n); i++) {
    177         if (radius->data.F32[i] > Rmax) {
     287    for (int i = 0; !done && (i < pixRadius->n); i++) {
     288        if (pixRadius->data.F32[i] > Rmax) {
    178289            // calculate the total flux for bin 'nOut'
    179290            float Area = M_PI*PS_SQR(Rmax);
    … …  
    185296                     nOut, radMax->data.F32[nOut], flux->data.F32[nOut], fluxErr->data.F32[nOut], fill->data.F32[nOut], Area);
    186297
     298            nPass ++;
     299            // if (nPass % 1000 == 0) {fprintf (stderr, "!");}
     300
    187301            nOut ++;
    188302            if (nOut >= radMax->n) break;
    … …  
    195309    flux->n = fluxErr->n = fill->n = nOut;
    196310   
    197     psFree(source->radial->flux);
    198     psFree(source->radial->fluxErr);
    199     psFree(source->radial->fill);
    200 
    201     source->radial->flux = flux;
    202     source->radial->fluxErr = fluxErr;
    203     source->radial->fill = fill;
    204 
    205     psFree (radius);
     311    radialAper->flux = flux;
     312    radialAper->fluxErr = fluxErr;
     313    radialAper->fill = fill;
     314
     315    psFree (pixRadius);
    206316    psFree (pixFlux);
    207317    psFree (pixVar);
    208318
     319    nCalls ++;
     320    // if (nCalls % 100 == 0) {fprintf (stderr, "*");}
    209321    return true;
    210322}

  • branches/eam_branches/ipp-20101205/psphot/src/psphotRadialAperturesByObject.c

    r29936 r30101  
    33// aperture-like measurements for extended sources
    44// flux in simple, circular apertures
    5 bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule) {
     5
     6// **** it looks like this function will re-point the source pixels at the specified FILERULE
     7// **** I need to distinguish PSF-matched images from raw
     8// **** save (somewhere) the PSF-matched PSF values
     9
     10// this function measures the radial aperture fluxes for the set of readouts.  this function
     11// may be called multiple times (presumably with different matched PSF sizes).  we must have
     12// already added an entry to the readout->analysis identifying the FWHM of this version.
     13
     14bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule, int nMatchedPSF) {
    615
    716    bool status;
    … …  
    2837    psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
    2938    psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
     39    float outerRadius = radMax->data.F32[radMax->n - 1];
    3040
    3141    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    3949    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
    4050
     51    // how many target PSFs do we want?
     52    int nPSFsizes = 0;
     53    {
     54        psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
     55        if (status) {
     56            nPSFsizes = 1;
     57        } else {
     58            psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
     59            psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
     60            nPSFsizes = fwhmValues->n;
     61        }
     62    }
     63   
    4164    // source analysis is done in S/N order (brightest first)
    4265    objects = psArraySort (objects, pmPhotObjSortBySN);
    … …  
    5376        psAssert (readout, "missing readout?");
    5477
     78        psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
     79        if (!fwhmValues) {
     80            psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
     81            return false;
     82        }
     83        if (fwhmValues->n != nMatchedPSF + 1) {
     84            psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM sequence is inconsistent this readout");
     85            return false;
     86        }
     87        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF FWHM of %s : %f pixels\n", file->name, fwhmValues->data.F32[nMatchedPSF]);
     88
    5589        readouts->data[i] = psMemIncrRefCounter(readout);
    5690    }
    5791
    5892    // process the objects in order. 
    59     for (int i = 0; i < objects->n; i++) {
     93    // XXX TEST: fix this!
     94    for (int i = 0; (i < 50) && (i < objects->n); i++) {
    6095        pmPhotObj *object = objects->data[i];
    6196        if (!object) continue;
    … …  
    77112            if (source->peak->SN < SN_LIM) continue;
    78113
     114            int index = source->imageID;
     115            if (index >= readouts->n) continue; // skip the sources generated by the chisq image
     116            pmReadout *readout = readouts->data[index];
     117
     118            // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
     119            // psphotVisualShowImage(readout);
     120
     121            // allocate pmSourceExtendedParameters, if not already defined
     122            if (!source->radialAper) {
     123                source->radialAper = psArrayAlloc(nPSFsizes);
     124            }
     125
    79126            // replace object in image
    80127            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    … …  
    83130            Nradial ++;
    84131
    85             int index = source->imageID;
    86             pmReadout *readout = readouts->data[index];
     132            // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
     133            // psphotVisualShowImage(readout);
    87134
    88135            // force source image to be a bit larger...
    89             float radius = source->peak->xf - source->pixels->col0;
    90             radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
    91             radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
    92             radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
    93             pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
     136            // float radius = source->peak->xf - source->pixels->col0;
     137            // radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
     138            // radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
     139            // radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
     140            pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
    94141
    95             if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) {
     142            if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, nMatchedPSF)) {
    96143                psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    97144            } else {
    … …  
    101148            // re-subtract the object, leave local sky
    102149            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     150
     151            // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
     152            // psphotVisualShowImage(readout);
    103153        }
    104154    }

  • branches/eam_branches/ipp-20101205/psphot/src/psphotSetThreads.c

    r29004 r30101  
    3535    psFree(task);
    3636
    37     task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 12);
     37    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 13);
    3838    task->function = &psphotExtendedSourceFits_Threaded;
    3939    psThreadTaskAdd(task);

  • branches/eam_branches/ipp-20101205/psphot/src/psphotStackChisqImage.c

    r30028 r30101  
    66
    77// XXX supply filename or keep PSPHOT.INPUT fixed?
    8 bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleCnv)
     8bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleSrc)
    99{
    1010    psTimerStart ("psphot.chisq.image");
    … …  
    2727
    2828    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.CHISQ.NUM", PS_META_REPLACE, "", num);
     29
     30    // we need to increment the counter for ruleDet and ruleSrc:
    2931    num++;
    30     psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", num);
    3132
    3233    // save the resulting image in the 'detection' set
    … …  
    3536        return false;
    3637    }
     38    psphotFileruleCountSet(config, ruleDet, num);
    3739
    38     // save the resulting image in the 'convolved' set
    39     // XXX whil am I doing this as well?
    40     if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleCnv, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
    41         psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files");
    42         return false;
     40    // also save the resulting image in the 'source' set (analysis set)
     41    if (strcmp(ruleDet, ruleSrc)) {
     42        if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleSrc, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
     43            psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files");
     44            return false;
     45        }
     46        psphotFileruleCountSet(config, ruleSrc, num);
    4347    }
    4448
    … …  
    120124    psAssert (status, "programming error: must define PSPHOT.CHISQ.NUM");
    121125
    122     int inputNum = psphotFileruleCount(config, "PSPHOT.INPUT");
     126    int inputNum = psphotFileruleCount(config, filerule);
    123127
    124128    pmFPAfileRemoveSingle (config->files, filerule, chisqNum);
    125129
    126130    inputNum --;
    127     psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", inputNum);
     131    psphotFileruleCountSet(config, filerule, inputNum);
    128132
    129133    return true;

  • branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFs.c

    r30006 r30101  
    106106    rescaleData(readoutOut, config, options, index);
    107107
    108     // dumpImage(readoutOut, readoutSrc, index, "convolved");
     108    // save the output fwhm values in the readout->analysis.  we may have / will have multiple output PSF sizes,
     109    // so we save this in a vector.  if the vector is not yet defined, create it
     110    bool mdok = false;
     111    psVector *fwhmValues = psMetadataLookupVector(&mdok, readoutOut->analysis, "STACK.PSF.FWHM.VALUES");
     112    if (!fwhmValues) {
     113        fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
     114        psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
     115        psFree(fwhmValues); // drops the extra copy
     116    }
     117    psVectorAppend(fwhmValues, options->targetSeeing);
    109118
    110119    return true;

  • branches/eam_branches/ipp-20101205/psphot/src/psphotStackPSF.c

    r28013 r30101  
    1212    bool autoPSF = psMetadataLookupBool (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.AUTO");
    1313
     14    // Get the recipe values
     15    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
     16    psAssert(recipe, "We've thrown an error on this before.");
     17       
     18    char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
     19
    1420    if (autoPSF) {
    15         // Get the recipe values
    16         psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
    17         psAssert(recipe, "We've thrown an error on this before.");
    18 
    1921        int psfInstances = psMetadataLookupS32(NULL, recipe, "PSF.INSTANCES"); // Number of instances for PSF
    2022        float psfRadius = psMetadataLookupF32(NULL, recipe, "PSF.RADIUS"); // Radius for PSF
    21         const char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
    2223        int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF
    2324
    … …  
    4950
    5051        float targetFWHM = psMetadataLookupF32 (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
    51         psAssert (isfinite(targetFWHM), "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
     52        if (!mdok) {
     53            psVector *fwhmValues = psMetadataLookupVector(&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
     54            psAssert (mdok, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
     55            targetFWHM = fwhmValues->data.F32[0];
     56        }
    5257
    5358        float Sxx = sqrt(2.0)*targetFWHM / 2.35;
    5459
    5560        // XXX probably should make the model type (and par 7) optional from recipe
    56         psf = pmPSFBuildSimple("PS_MODEL_PS1_V1", Sxx, Sxx, 0.0, 1.0);
     61        psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
    5762        if (!psf) {
    5863            psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF.");

  • branches/eam_branches/ipp-20101205/psphot/src/psphotStackParseCamera.c

    r29548 r30101  
    1515    }
    1616
     17    int nRaw = 0;
     18    int nCnv = 0;
    1719    int nInputs = inputs->list->n;
    1820    for (int i = 0; i < nInputs; i++) {
    … …  
    5759                }
    5860            }
     61            nRaw ++;
    5962        }
    6063
    … …  
    8891                }
    8992            }
     93            nCnv ++;
    9094        }
    9195
    … …  
    9397            psError(PSPHOT_ERR_CONFIG, true, "Component %s (%d) lacks both RAW:IMAGE and CNV:IMAGE of type STR", item->name, i);
    9498            return false;
     99        }
     100
     101        // XXX what if they do not match in length
     102        if (nCnv && nRaw) {
     103            if (nCnv != nRaw) {
     104                psError (PSPHOT_ERR_CONFIG, true, "if both RAW and CNV images are supplied, the number must match");
     105                return false;
     106            }
    95107        }
    96108
    … …  
    150162    }
    151163    psMetadataRemoveKey(config->arguments, "FILENAMES");
     164    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.INPUT.RAW.NUM", PS_META_REPLACE, "number of inputs", nRaw);
     165    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.INPUT.CNV.NUM", PS_META_REPLACE, "number of inputs", nCnv);
     166    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.OUTPUT.IMAGE.NUM", PS_META_REPLACE, "number of inputs", nInputs);
    152167    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", nInputs);
    153168

  • branches/eam_branches/ipp-20101205/psphot/src/psphotStackReadout.c

    r30027 r30101  
    11# include "psphotInternal.h"
    22
     3// we have 3 possible real filesets:
    34# define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
    45# define STACK_CNV "PSPHOT.STACK.INPUT.CNV"
    5 # define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"
     6# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"  /* the psf-matched image */
     7
     8// we have 3 files on which we operate:
     9// DET (detection image)       : nominally RAW (optionally CNV?)
     10// SRC (source analysis image) : nominally CNV (optionally RAW)
     11// OUT (psf-matched images)    : always OUT
    612
    713bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
    … …  
    2127    PS_ASSERT_PTR_NON_NULL (breakPt, false);
    2228
    23     BST = RAW : CNV;
     29    // XXX or do I set OUT to be a pmFPAfile pointing at the input of interest?
     30    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
     31    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
     32    char *STACK_DET = STACK_RAW; // XXX optionally allow this to be CNV?
    2433
    2534    // we have 3 relevant files: RAW, CNV, OUT
    2635
    2736    // set the photcode for each image
    28     if (!psphotAddPhotcode (config, view, STACK_BST)) {
     37    if (!psphotAddPhotcode (config, view, STACK_SRC)) {
    2938        psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
    3039        return false;
    … …  
    3342    // Generate the mask and weight images
    3443    // XXX this should be done before we perform the convolutions
    35     if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) {
    36         return psphotReadoutCleanup (config, view, STACK_OUT);
     44    if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
     45        return psphotReadoutCleanup (config, view, STACK_SRC);
    3746    }
    3847    if (!strcasecmp (breakPt, "NOTHING")) {
    39         return psphotReadoutCleanup (config, view, STACK_OUT);
     48        return psphotReadoutCleanup (config, view, STACK_SRC);
    4049    }
    4150
    … …  
    4352    // XXX I think this is not defined correctly for an array of images.
    4453    // XXX probably need to subtract the model (same model?) for both RAW and OUT
    45     if (!psphotModelBackground (config, view, STACK_RAW)) {
    46         return psphotReadoutCleanup (config, view, STACK_OUT);
    47     }
    48     if (!psphotSubtractBackground (config, view, STACK_RAW)) {
    49         return psphotReadoutCleanup (config, view, STACK_OUT);
     54    if (!psphotModelBackground (config, view, STACK_DET)) {
     55        return psphotReadoutCleanup (config, view, STACK_SRC);
     56    }
     57    if (!psphotSubtractBackground (config, view, STACK_DET)) {
     58        return psphotReadoutCleanup (config, view, STACK_SRC);
    5059    }
    5160    if (!strcasecmp (breakPt, "BACKMDL")) {
    52         return psphotReadoutCleanup (config, view, STACK_OUT);
    53     }
    54 
    55     // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are determined and saved on
    56     // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT
    57     if (!psphotLoadPSF (config, view, STACK_RAW)) {
    58         psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
    59         return psphotReadoutCleanup (config, view, STACK_OUT);
    60     }
    61 
    62     if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_OUT)) {
     61        return psphotReadoutCleanup (config, view, STACK_SRC);
     62    }
     63
     64    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
    6365        psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
    64         return psphotReadoutCleanup (config, view, STACK_OUT);
     66        return psphotReadoutCleanup (config, view, STACK_SRC);
    6567    }
    6668    if (!strcasecmp (breakPt, "CHISQ")) {
    67         return psphotReadoutCleanup (config, view, STACK_OUT);
     69        return psphotReadoutCleanup (config, view, STACK_SRC);
    6870    }
    6971
    7072    // find the detections (by peak and/or footprint) in the image.
    7173    // This finds the detections on Chisq image as well as the individuals
    72     if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
     74    if (!psphotFindDetections (config, view, STACK_DET, true)) { // pass 1
    7375        // this only happens if we had an error in psphotFindDetections
    7476        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    75         return psphotReadoutCleanup (config, view, STACK_OUT);
    76     }
    77 
    78     // copy the detections from RAW to OUT
    79     if (!psphotCopySources (config, view, STACK_OUT, STACK_RAW)) {
    80         psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    81         return psphotReadoutCleanup (config, view, STACK_OUT);
     77        return psphotReadoutCleanup (config, view, STACK_SRC);
     78    }
     79
     80    // copy the detections from DET to SRC
     81    if (strcmp(STACK_SRC, STACK_DET)) {
     82        if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
     83            psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
     84            return psphotReadoutCleanup (config, view, STACK_SRC);
     85        }
    8286    }
    8387
    8488    // construct sources and measure basic stats (saved on detections->newSources)
    85     // only run this on detections from the input images, not chisq image
    86     if (!psphotSourceStats (config, view, STACK_OUT, true)) { // pass 1
     89    if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
    8790        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    88         return psphotReadoutCleanup (config, view, STACK_OUT);
     91        return psphotReadoutCleanup (config, view, STACK_SRC);
    8992    }
    9093
    9194    // generate the objects (object unify the sources from the different images)
    92     psArray *objects = psphotMatchSources (config, view, STACK_OUT);
     95    // XXX this could just match the detections for the chisq image, and not bother measuring the
     96    // source stats in that case...
     97    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
    9398
    9499    // construct sources for the newly-generated sources (from other images)
    95     if (!psphotSourceStats (config, view, STACK_OUT, false)) { // pass 1
     100    if (!psphotSourceStats (config, view, STACK_SRC, false)) { // pass 1
    96101        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    97         return psphotReadoutCleanup (config, view, STACK_OUT);
     102        return psphotReadoutCleanup (config, view, STACK_SRC);
    98103    }
    99104
    … …  
    101106    // if (!psphotDeblendSatstars (config, view)) {
    102107    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
    103     //     return psphotReadoutCleanup (config, view, STACK_OUT);
     108    //     return psphotReadoutCleanup (config, view, STACK_SRC);
    104109    // }
    105110
    … …  
    107112    // if (!psphotBasicDeblend (config, view)) {
    108113    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
    109     //     return psphotReadoutCleanup (config, view, STACK_OUT);
     114    //     return psphotReadoutCleanup (config, view, STACK_SRC);
    110115    // }
    111116
    112117    // classify sources based on moments, brightness
    113118    // only run this on detections from the input images, not chisq image
    114     if (!psphotRoughClass (config, view, STACK_OUT)) {
     119    if (!psphotRoughClass (config, view, STACK_SRC)) {
    115120        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
    116         return psphotReadoutCleanup (config, view, STACK_OUT);
     121        return psphotReadoutCleanup (config, view, STACK_SRC);
    117122    }
    118123    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
    119124    // only run this on detections from the input images, not chisq image
    120     if (!psphotImageQuality (config, view, STACK_OUT)) { // pass 1
     125    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
    121126        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
    122         return psphotReadoutCleanup (config, view, STACK_OUT);
     127        return psphotReadoutCleanup (config, view, STACK_SRC);
    123128    }
    124129    if (!strcasecmp (breakPt, "MOMENTS")) {
    125         return psphotReadoutCleanup (config, view, STACK_OUT);
     130        return psphotReadoutCleanup (config, view, STACK_SRC);
    126131    }
    127132
    128133    // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
    129134    // this step is skipped
    130     if (!psphotChoosePSF (config, view, STACK_OUT)) { // pass 1
     135    if (!psphotChoosePSF (config, view, STACK_SRC)) { // pass 1
    131136        psLogMsg ("psphot", 3, "failure to construct a psf model");
    132         return psphotReadoutCleanup (config, view, STACK_OUT);
     137        return psphotReadoutCleanup (config, view, STACK_SRC);
    133138    }
    134139    if (!strcasecmp (breakPt, "PSFMODEL")) {
    135         return psphotReadoutCleanup (config, view, STACK_OUT);
     140        return psphotReadoutCleanup (config, view, STACK_SRC);
    136141    }
    137142
    138143    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    139     psphotGuessModels (config, view, STACK_OUT);
     144    psphotGuessModels (config, view, STACK_SRC);
    140145
    141146    // merge the newly selected sources into the existing list
    142147    // NOTE: merge OLD and NEW
    143     psphotMergeSources (config, view, STACK_OUT);
     148    psphotMergeSources (config, view, STACK_SRC);
    144149
    145150    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    … …  
    147152
    148153    // identify CRs and extended sources
    149     psphotSourceSize (config, view, STACK_OUT, TRUE);
     154    psphotSourceSize (config, view, STACK_SRC, TRUE);
    150155
    151156    // measure aperture photometry corrections
    152     if (!psphotApResid (config, view, STACK_OUT)) {
     157    if (!psphotApResid (config, view, STACK_SRC)) {
    153158        psFree (objects);
    154159        psLogMsg ("psphot", 3, "failed on psphotApResid");
    155         return psphotReadoutCleanup (config, view, STACK_OUT);
     160        return psphotReadoutCleanup (config, view, STACK_SRC);
    156161    }
    157162
    158163    psphotStackObjectsUnifyPosition (objects);
    159164
    160     // measure circular, radial apertures (objects sorted by S/N)
    161     psphotRadialAperturesByObject (config, objects, view, STACK_OUT);
    162 
    163165    // measure elliptical apertures, petrosians (objects sorted by S/N)
    164     psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_OUT); // pass 1 (detections->allSources)
     166    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
    165167
    166168    // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
    167     psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources)
     169    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
    168170
    169171    // calculate source magnitudes
    170     psphotMagnitudes(config, view, STACK_OUT);
     172    psphotMagnitudes(config, view, STACK_SRC);
     173
     174    // copy the detections from SRC to OUT (for radial aperture photometry)
     175    if (!psphotCopySources (config, view, STACK_OUT, STACK_SRC)) {
     176        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
     177        return psphotReadoutCleanup (config, view, STACK_SRC);
     178    }
     179    // XXX need to do something to reassign the source pixels here
     180
     181    bool smoothAgain = true;
     182    for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
     183        // measure circular, radial apertures (objects sorted by S/N)
     184        psphotRadialAperturesByObject (config, objects, view, STACK_OUT, nMatchedPSF);
     185        psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
     186    }
    171187
    172188    if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
    … …  
    179195
    180196    // replace background in residual image
    181     psphotSkyReplace (config, view, STACK_RAW);
     197    psphotSkyReplace (config, view, STACK_DET);
    182198
    183199    // drop the references to the image pixels held by each source
    184     psphotSourceFreePixels (config, view, STACK_OUT);
    185 
    186     // remove chisq image from config->file:PSPHOT.INPUT (why?)
    187     psphotStackRemoveChisqFromInputs(config, STACK_RAW);
     200    // psphotSourceFreePixels (config, view, STACK_OUT);
     201    psphotSourceFreePixels (config, view, STACK_SRC);
     202
     203    // remove chisq image from config->file:PSPHOT.INPUT
     204    psphotStackRemoveChisqFromInputs(config, STACK_DET);
     205    if (strcmp(STACK_SRC, STACK_DET)) {
     206        psphotStackRemoveChisqFromInputs(config, STACK_SRC);
     207    }
    188208
    189209    psFree (objects);
    190210
    191211    // create the exported-metadata and free local data
    192     return psphotReadoutCleanup (config, view, STACK_OUT);
     212    return psphotReadoutCleanup (config, view, STACK_SRC);
    193213}
    194214
    … …  
    258278 ******
    259279
    260  the above is all wrong:  first, we sohould be doing the full
     280 the above is all wrong:  first, we should be doing the full
    261281 morphology analysis (ExtendedAnalysis & ExtendedFits) on the CNV or
    262282 RAW image (as desired optionally), etc.
Note: See TracChangeset for help on using the changeset viewer.