IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 31158 for trunk/ppStack


Ignore:
Timestamp:
Apr 4, 2011, 1:16:41 PM (15 years ago)
Author:
eugene
Message:

improvements to the output headers; some code organization; plug leaks; user-specified limits to input seeing ranges

Location:
trunk/ppStack
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack

    • Property svn:ignore

      •  

        old new  
        1717autom4te.cache
        1818Doxyfile
         19a.out.dSYM
    • Property svn:mergeinfo deleted

  • trunk/ppStack/src/ppStackArguments.c

    r30620 r31158  
    110110{
    111111    assert(config);
     112
     113    // generic arguments (version, dumpconfig)
     114    PS_ARGUMENTS_GENERIC( ppStack, config, argc, argv );
     115
     116    // thread arguments
     117    PS_ARGUMENTS_THREADS( ppStack, config, argc, argv )
    112118
    113119    {
    … …  
    181187    psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-variance", 0, "Suffix for temporary variance maps", NULL);
    182188    psMetadataAddBool(arguments, PS_LIST_TAIL, "-temp-delete", 0, "Delete temporary files on completion?", false);
    183     psMetadataAddS32(arguments, PS_LIST_TAIL, "-threads", 0, "Number of threads to use", 0);
    184189    psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "visualisation", false);
    185190
    … …  
    238243    valueArgStr(arguments, "-stats", "STATS", arguments);
    239244
    240     int numThreads = psMetadataLookupS32(NULL, arguments, "-threads"); // Number of threads
    241     if (numThreads > 0 && !psThreadPoolInit(numThreads)) {
    242         psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to setup %d threads", numThreads);
    243         return false;
    244     }
    245 
    246245    psMetadataAddBool(arguments, PS_LIST_TAIL, "PPSTACK.DEBUG.STACK", 0,
    247246                      "Read old convolved images to debug stack?", debugStack);

  • trunk/ppStack/src/ppStackCombineFinal.c

    r31054 r31158  
    6161        }
    6262
    63         // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
     63        // calls ppStackReadoutFinal(config, outRO, readouts, rejected) in ppStackReadout.c
    6464        psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
    6565        psArrayAdd(job->args, 1, thread);

  • trunk/ppStack/src/ppStackConvolve.c

    r31054 r31158  
    44
    55// Update the value of a concept
    6 #define UPDATE_CONCEPT(SOURCE, NAME, VALUE) { \
    7     psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \
    8     psAssert(item, "Concept should be present"); \
    9     psAssert(item->type == PS_TYPE_F32, "Concept should be F32"); \
    10     item->data.F32 = VALUE; \
    11 }
     6#define UPDATE_CONCEPT(SOURCE, NAME, VALUE) {                           \
     7        psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \
     8        psAssert(item, "Concept should be present");                    \
     9        psAssert(item->type == PS_TYPE_F32, "Concept should be F32");   \
     10        item->data.F32 = VALUE;                                         \
     11    }
    1212
    1313bool ppStackConvolve(ppStackOptions *options, pmConfig *config)
    … …  
    4747    psVectorInit(renorms, NAN);
    4848
     49    psVector *satValues = psVectorAllocEmpty(num, PS_TYPE_F32); // Renormalisation values for variances
     50
    4951    psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
    5052    psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
    … …  
    7072            psFree(cellList);
    7173            psFree(target);
     74            psFree(renorms);
     75            psFree(satValues);
    7276            return false;
    7377        }
    … …  
    8791            psFree(cellList);
    8892            psFree(target);
     93            psFree(renorms);
     94            psFree(satValues);
    8995            return false;
    9096        }
    … …  
    106112                psFree(cellList);
    107113                psFree(target);
     114                psFree(renorms);
     115                psFree(satValues);
    108116                return false;
    109117                // Non-fatal errors
    … …  
    120128                    psFree(cellList);
    121129                    psFree(target);
     130                    psFree(renorms);
     131                    psFree(satValues);
    122132                    return false;
    123133                }
    … …  
    166176            psFree(rng);
    167177            psFree(target);
     178            psFree(renorms);
     179            psFree(satValues);
    168180            return false;
    169181        }
    … …  
    177189            psFree(maskHeader);
    178190            psFree(target);
     191            psFree(renorms);
     192            psFree(satValues);
    179193            return false;
    180194        }
    … …  
    186200            psFree(rng);
    187201            psFree(target);
     202            psFree(renorms);
     203            psFree(satValues);
    188204            return false;
    189205        }
    … …  
    221237        if (options->matchZPs) {
    222238            // I think I need to take off the exposure time because we're going to set the new exposure time
     239            // Clarification: the zero point (ZP) in the header should be set such that:
     240            // M_app = m_inst + ZP + 2.5*log(exptime), where exptime in the output is sumExposure
    223241            psMetadataItem *zpItem = psMetadataLookup(inCell->parent->parent->concepts, "FPA.ZP");
     242            float inZP = zpItem->data.F32;
    224243            zpItem->data.F32 += options->norm->data.F32[i] + 2.5*log10(options->sumExposure);
    225244
    … …  
    228247
    229248            expItem = psMetadataLookup(inCell->concepts, "CELL.EXPOSURE");
     249            float inExptime = expItem->data.F32;
    230250            expItem->data.F32 = options->sumExposure;
     251
     252            // flux_out = flux_in * ten(-0.4*norm) -- save the individual saturation values
     253            psMetadataItem *satItem = psMetadataLookup(inCell->concepts, "CELL.SATURATION");
     254            float inSat = satItem->data.F32;
     255            satItem->data.F32 *= pow(10.0, -0.4*options->norm->data.F32[i]);
     256            psVectorAppend (satValues, satItem->data.F32);
     257
     258            psLogMsg("ppStack", PS_LOG_INFO, "image %d mods : zp %f -> %f, exptime %f -> %f, sat %f -> %f",
     259                     i, inZP, zpItem->data.F32, inExptime, expItem->data.F32, inSat, satItem->data.F32);
    231260        }
    232261
    … …  
    237266            psFree(rng);
    238267            psFree(target);
     268            psFree(renorms);
     269            psFree(satValues);
    239270            return false;
    240271        }
    … …  
    257288        psFree(fpaList);
    258289        psFree(cellList);
     290        psFree(renorms);
     291        psFree(satValues);
    259292        return true;
    260293    }
    … …  
    281314        psFree(fpaList);
    282315        psFree(cellList);
    283 
    284         UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",  options->sumExposure);
    285         UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", options->sumExposure);
    286         UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
    287         UPDATE_CONCEPT(outFPA,  "FPA.ZP",        options->zp);
    288 
    289         UPDATE_CONCEPT(unconvFPA,  "FPA.EXPOSURE",  options->sumExposure);
    290         UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE", options->sumExposure);
    291         UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME", NAN);
    292         UPDATE_CONCEPT(unconvFPA,  "FPA.ZP",        options->zp);
     316       
     317        // The best guess for an output saturation value depends on the recipe.  If we have
     318        // 'safe' on, the we require at least 2 pixels to generate a valid output pixel.  In
     319        // this case, the best value for CELL.SATURATION is the 2nd highest value in the list.
     320        // If not, it should be the higest value in the list
     321        bool mdok = false;
     322        bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
     323        psVectorSortInPlace(satValues);
     324        float satBest = safe && satValues->n > 1 ? satValues->data.F32[1] : satValues->data.F32[0];
     325
     326        // UPDATE CELL.SATURATION here
     327        UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",    options->sumExposure);
     328        UPDATE_CONCEPT(outCell, "CELL.EXPOSURE",   options->sumExposure);
     329        UPDATE_CONCEPT(outCell, "CELL.DARKTIME",   NAN);
     330        UPDATE_CONCEPT(outCell, "CELL.SATURATION", satBest);
     331        UPDATE_CONCEPT(outFPA,  "FPA.ZP",          options->zp);
     332        UPDATE_CONCEPT(outFPA,  "FPA.AIRMASS",     options->airmass);
     333
     334        UPDATE_CONCEPT(unconvFPA,  "FPA.EXPOSURE",    options->sumExposure);
     335        UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE",   options->sumExposure);
     336        UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME",   NAN);
     337        UPDATE_CONCEPT(unconvCell, "CELL.SATURATION", satBest);
     338        UPDATE_CONCEPT(unconvFPA,  "FPA.ZP",          options->zp);
     339        UPDATE_CONCEPT(unconvFPA,  "FPA.AIRMASS",     options->airmass);
     340
     341        psLogMsg("ppStack", PS_LOG_INFO, "stack adjust metadata values : zp %f, exptime %f, sat %f", options->zp, options->sumExposure, satBest);
    293342
    294343        if (options->stats) {
    … …  
    297346            double time = psTimeToMJD(fpaTime);
    298347            psMetadataAddF64(options->stats, PS_LIST_TAIL, "MJD_OBS", PS_META_DUPLICATE_OK,
    299                                  "Average MJD_OBS of inputs", time);
    300 
    301         }
    302     }
    303 
     348                             "Average MJD_OBS of inputs", time);
     349
     350        }
     351    }
     352
     353    // XXX EAM : this may be overly harsh -- or at least it would be if I (EAM) hadn't changed
     354    // the values of matchChi2 I modified pmSubtraction.c to fit a 2nd order polynomial to the
     355    // star chisq distribution (because of systematic errors in the model being matched).  This
     356    // fit forces the mean value to be 0.0.  Perhaps we can / should exclude images which have
     357    // an excessively high value for the rms?
    304358
    305359    // Reject images out-of-hand on the basis of their match chi^2
    … …  
    316370            psError(PPSTACK_ERR_PROG, false, "Unable to sort vector.");
    317371            psFree(values);
     372            psFree(renorms);
     373            psFree(satValues);
    318374            return false;
    319375        }
    … …  
    356412        psErrorClear();
    357413        psWarning("No good images survived convolution stage.");
     414        psFree(renorms);
     415        psFree(satValues);
    358416        return true;
    359417    }
    … …  
    366424    }
    367425    psFree(renorms);
     426    psFree(satValues);
    368427
    369428    if (options->stats) {

  • trunk/ppStack/src/ppStackFiles.c

    r30620 r31158  
    6565    for (int i = 0; i < numLeaks; i++) {
    6666        psMemBlock *mb = leaks[i];
    67         fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize,
    68                 mb->file, mb->lineno);
     67        fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize, mb->file, mb->lineno);
    6968        total += mb->userMemorySize;
    7069    }
    … …  
    7574    num++;
    7675}
    77 
    78 
    7976
    8077// Activate/deactivate a list of files

  • trunk/ppStack/src/ppStackLoop.c

    r30722 r31158  
    5555    }
    5656
    57     // Start threading
    58     ppStackThreadInit();
     57    // Define threading elements
    5958    ppStackThreadData *stack = ppStackThreadDataSetup(options, config, true);
    6059    if (!stack) {
    … …  
    114113    }
    115114
    116     // Final combination
     115    // Final combination.  This one does NOT need to be normalized since the convolution takes care of that
    117116    psTrace("ppStack", 2, "Final stack of convolved images....\n");
    118117    if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, false, true)) {
    … …  
    181180        }
    182181
    183         // generate the unconvolved stack
     182        // generate the unconvolved stack. NOTE: this one must be normalized since the inputs have not been
    184183        psTrace("ppStack", 2, "Stack of unconvolved images....\n");
    185         if (!ppStackCombineFinal(stack, options->origCovars, options, config,
    186                                  false, true, false)) {
     184        if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false)) {
    187185            psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination.");
    188186            psFree(stack);

  • trunk/ppStack/src/ppStackMatch.c

    r30685 r31158  
    142142    bool mdok;                          // Status of MD lookup
    143143    float penalty = psMetadataLookupF32(NULL, ppsub, "PENALTY"); // Penalty for wideness
    144     int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
     144    int threads = psMetadataLookupS32(NULL, config->arguments, "NTHREADS"); // Number of threads
    145145
    146146    // Replaced pmReadoutMaskNonfinite with pmReadoutMaskInvalid (tests for already masked pixels)

  • trunk/ppStack/src/ppStackOptions.h

    r30620 r31158  
    2222    float sumExposure;                  // Sum of exposure times
    2323    float zp;                           // Zero point for output
     24    float zpErr;                        // Zero point error for output
     25    float airmass;                      // airmass for output image
     26    float airmassSlope;                 // airmass slope used for zp analysis
    2427    psVector *inputMask;                // Mask for inputs
    2528    psArray *sourceLists;               // Individual lists of sources for matching
    2629    psVector *norm;                     // Normalisation for each image
     30    psVector *zpInput;                  // reported zero point of input exposure
     31    psVector *expTimeInput;             // exposure time of input exposure
     32    psVector *airmassInput;             // airmass of input exposure
    2733    psArray *sources;                   // Matched sources
    2834    // Convolve

  • trunk/ppStack/src/ppStackPrepare.c

    r30620 r31158  
    1 #include "ppStack.h"
    2 
    3 #define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
    4                           PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
    5                           PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources
     1# include "ppStack.h"
     2
     3# define RE_PHOTOMETER 0
     4
     5# define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
     6                           PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
     7                           PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources
    68
    79bool ppStackInputPhotometer(const pmReadout *ro, const psArray *sources, const pmConfig *config)
    … …  
    1315    psAssert(recipe, "We've thrown an error on this before.");
    1416
     17# if (RE_PHOTOMETER)
    1518    bool mdok;                          // Status of MD lookup
    1619
    … …  
    3841    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask
    3942
     43    // Measure background
    4044    psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout
    41 
    42     // Measure background
    4345    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    4446    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
    … …  
    5254    psFree(stats);
    5355    psFree(rng);
    54 
    55     // Photometer sources
    56     int numCols = image->numCols, numRows = image->numRows; // Size of image
     56# endif
     57
     58    // Examine the sources (photometer if RE_PHOTOMETER is set)
    5759    int numSources = sources->n; // Number of sources
    5860    int numGood = 0;            // Number of good sources
    … …  
    6466            continue;
    6567        }
    66         float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
     68
     69// This code re-measures the magnitude in the assumption that the warp analysis failed It would
     70// be better to catch and correct such cases.  In any case, we should check on the validity of
     71// this assumption (eg, compare the mean psf-ap offset) and raise an error on that?
     72
     73# if (RE_PHOTOMETER)
     74       
     75        int numCols = image->numCols, numRows = image->numRows; // Size of image
     76        float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
    6777        int xCentral = x + 0.5, yCentral = y + 0.5; // Central pixel
    6878        // Range of integration
    … …  
    96106            }
    97107        }
     108# else
     109        numGood++;
     110        maxMag = PS_MAX(maxMag, source->psfMag);
     111# endif
    98112    }
    99113    psLogMsg("ppStack", PS_LOG_INFO, "Photometered %d good sources in image; max mag %f", numGood, maxMag);
    … …  
    217231        }
    218232    }
     233
     234    bool mdok = false;
     235    float maxFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.MAX"); // max allowed input fwhm
     236    float clipFWHMnSig = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.NSIGMA"); // sigma clipping of inputs
    219237
    220238    psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
    … …  
    251269        }
    252270
     271        // reject any input images which exceed the specified max FWHM
     272        if (isfinite(maxFWHM) && (options->inputSeeing->data.F32[i] > maxFWHM)) {
     273            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
     274            psLogMsg("ppStack", PS_LOG_INFO, "PSF FWHM for image %d is too large (%f vs %f maxFWHM) --- rejected.", i, options->inputSeeing->data.F32[i], maxFWHM);
     275        }
     276
    253277        psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]);
    254278    }
    … …  
    257281    }
    258282    psFree(log);
     283
     284    // We should have the ability to filter the input list based on the seeing:
     285    // * reject above some max value and/or min value
     286    // * reject images out-of-line with the rest of the images
     287
     288    // measure stats
     289    psStats *fwhmStats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     290    psVectorStats (fwhmStats, options->inputSeeing, NULL, options->inputMask, 0xff);
     291    psLogMsg("ppStack", PS_LOG_INFO, "Input FWHMs : %f +/- %f", fwhmStats->clippedMean, fwhmStats->clippedStdev);
     292
     293    if (isfinite(clipFWHMnSig)) {
     294        float fwhmLimit = fwhmStats->clippedMean + clipFWHMnSig * fwhmStats->clippedStdev;
     295        fwhmLimit = isfinite(maxFWHM) ? PS_MIN (maxFWHM, fwhmLimit) : fwhmLimit;
     296        for (int i = 0; i < num; i++) {
     297            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
     298            if (options->inputSeeing->data.F32[i] > fwhmLimit) {
     299                options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
     300                psLogMsg("ppStack", PS_LOG_INFO, "PSF FWHM for image %d is too large (%f vs %f fwhmLimit) --- rejected.", i, options->inputSeeing->data.F32[i], fwhmLimit);
     301            }
     302        }
     303    }
    259304
    260305    // Generate target PSF

  • trunk/ppStack/src/ppStackReject.c

    r30620 r31158  
    5555
    5656    // Reject bad pixels
    57     if (psMetadataLookupS32(NULL, config->arguments, "-threads") > 0) {
     57    if (psMetadataLookupS32(NULL, config->arguments, "NTHREADS") > 0) {
    5858        pmStackRejectThreadsInit();
    5959    }

  • trunk/ppStack/src/ppStackSetup.c

    r30620 r31158  
    5151    psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments,
    5252                                                           "OUTPUT")); // Name for temporary files
    53     const char *tempName = basename(outputName);
     53    const char *tempName = psStringFileBasename(outputName);
    5454    if (!tempName) {
    5555        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to construct basename for temporary files.");
    … …  
    6565                "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
    6666        psFree(outputName);
     67        psFree(tempName);
    6768        return false;
    6869    }
    … …  
    8283    }
    8384    psFree(outputName);
     85    psFree(tempName);
    8486
    8587    // Original images

  • trunk/ppStack/src/ppStackSources.c

    r30685 r31158  
    115115    float fracMatch = psMetadataLookupF32(NULL, recipe, "ZP.MATCH"); // Fraction of images to match for star
    116116
    117     psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms
     117    bool mdok = false;
     118    float airmassTarget = psMetadataLookupF32(&mdok, recipe, "ZP.AIRMASS.TARGET"); // output airmass value
     119    if (!mdok) {
     120        airmassTarget = 1.0;
     121    }
     122
     123    psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms (slopes) by filter
    118124    if (!airmassZP) {
    119125        psError(PPSTACK_ERR_CONFIG, false, "Unable to find ZP.AIRMASS in recipe.");
    … …  
    128134    int num = options->num;             // Number of inputs
    129135    psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n);
     136
     137    // vectors to store these inputs so they may be recorded in the output headers
     138    options->zpInput      = psVectorAlloc(num, PS_TYPE_F32);
     139    options->expTimeInput = psVectorAlloc(num, PS_TYPE_F32);
     140    options->airmassInput = psVectorAlloc(num, PS_TYPE_F32);
    130141
    131142    psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Relative zero points for each image
    … …  
    167178        const char *expFilter = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.FILTER"); // Filter name
    168179        zpExp->data.F32[i] = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.ZP"); // Exposure zero point
     180        // XXX need to get the zero point error values and propagate to get the FPA.ZP.ERR value
     181
     182        options->zpInput->data.F32[i] = zpExp->data.F32[i]; // NOTE zpExp may be re-assigned below using relative photometry
     183        options->expTimeInput->data.F32[i] = exptime;
     184        options->airmassInput->data.F32[i] = airmass;
     185
    169186        psLogMsg("ppStack", PS_LOG_INFO,
    170187                 "Image %d: %.2f sec exposure in %s at airmass %.2f with zero point %.2f",
    … …  
    185202        if (!filter) {
    186203            filter = expFilter;
    187             bool mdok;
    188204            airmassTerm = psMetadataLookupF32(&mdok, airmassZP, filter);
    189205            if (!mdok || !isfinite(airmassTerm)) {
    … …  
    206222        }
    207223
     224        // XXX this is wrong, or at least inconsistent with the above: this needs to include
     225        // a value for the nominal system zero point to be consistent with zpExp
    208226        zp->data.F32[i] = airmassTerm * airmass + 2.5 * log10(exptime);
    209227    }
    … …  
    236254
    237255    if (zpExpNum == numGoodImages) {
     256        psLogMsg ("ppStack", PS_LOG_INFO, "all zero points are finite; using reported zero points listed above");
    238257        for (int i = 0; i < num; i++) {
    239258            zp->data.F32[i] = zpExp->data.F32[i];
     259        }
     260    } else {
     261        psLogMsg ("ppStack", PS_LOG_INFO, "missing some zero points; using guess values:");
     262        for (int i = 0; i < num; i++) {
     263            psLogMsg("ppStack", PS_LOG_INFO, "Image %d: %.2f sec exposure with zero point %.2f", i, options->exposures->data.F32[i], zp->data.F32[i]);
    240264        }
    241265    }
    … …  
    265289            }
    266290        }
    267         // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
    268         // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
    269         // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
    270         // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
    271         // We don't need to know the magnitude zero point for the filter, since it cancels out
     291
     292        // EAM : the discussion here was not quite right (or at least sloppy).  Here is a replacement explanation:
     293
     294        // For any star, the observed instrumental magnitude on an image and the apparent magnitude are related by:
     295        // M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency
     296        // NOTE the sign of 'transparency'  this must agree with the definition in pmSourceMatch.c. see, eg, line 457 where
     297        // transparency = m_inst + zp + c1 * airmass + 2.5log(t) - M_app
     298
     299        // we want to adjust the input images to be in a consistent flux system so that the
     300        // final stack can be generated with a specific target zero point.  Any adjustment to
     301        // the flux scale of the image must be made in coordination with the resulting
     302        // zeropoint, exposure time, and airmass such that the above relationship yields the
     303        // same apparent magnitude for a given star:
     304
     305        // m_inst_i : instrumental mags on input image (in)
     306        // m_inst_o : instrumental mags on re-normalized image (out)
     307
     308        // m_inst_o + zp_o + c1 * airmass_o + 2.5log(t_o) - trans_o = m_inst_i + zp_i + c1 * airmass_i + 2.5log(t_i) - trans_i
     309
     310        // m_inst_o = m_inst_i + (zp_i - zp_o) + c1 * (airmass_i - airmass_o) + 2.5log(t_i) - 2.5log(t_o) - trans_i + trans_o
     311
     312        // zp_i, airmass_i, t_i, trans_i : reported or measured for input image
     313
     314        // zp_o      = zpTarget      (from recipe)
     315        // airmass_o = airmassTarget (from recipe)
     316        // t_o       = sumExpTime    [sum of input exposure times: once images are scale to this time, they can be avereaged]
     317        // trans_o   = 0.0           [obviously!]
     318
     319        // we have 2 cases: (a) all reported ZPs are good or (b) some are bad:
     320        // (a) FPA.ZP = zp_i + c1 * airmass_i
     321        //  --> zp[i] = zp_i + c1 * airmass_i + 2.5log(exptime_i)
     322        // (b)  zp[i] = c1 * airmass_i + 2.5log(exptime_i)
     323        // NOTE: in case (b), the current code is equating the TARGET zp with the NOMINAL zp, which is wrong.
     324
     325        // m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i
     326
    272327        if (options->matchZPs) {
    273328            options->norm = psVectorAlloc(num, PS_TYPE_F32);
    … …  
    277332                }
    278333                psArray *sources = sourceLists->data[i]; // Sources of interest
    279                 float magCorr = zp->data.F32[i] - trans->data.F32[i] - 2.5*log10(options->sumExposure);
    280                 if (zpExpNum == numGoodImages) {
     334                float magCorr = zp->data.F32[i] - trans->data.F32[i] - 2.5*log10(options->sumExposure) - airmassTerm * airmassTarget;
     335                if (zpExpNum == numGoodImages) { // case (a)
    281336                    // Using measured zero points, so attempt to set target zero point
     337                    // XXX see NOTE above regarding case (b) : this is wrong.  the code should load a nominal zero point and supply it above
     338                    //
    282339                    magCorr -= zpTarget;
    283340                }
    … …  
    302359            // Producing image with target zero point
    303360            options->zp = zpTarget;
     361            options->airmass = airmassTarget;
     362            options->airmassSlope = airmassTerm;
    304363        } else {
    305364            options->zp = NAN;

  • trunk/ppStack/src/ppStackThread.c

    r30620 r31158  
    101101    }
    102102
    103     int numThreads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
     103    int numThreads = psMetadataLookupS32(NULL, config->arguments, "NTHREADS"); // Number of threads
    104104
    105105    // Generate readouts for each input file in each file group
    … …  
    248248
    249249
    250 void ppStackThreadInit(void)
     250void ppStackSetThreads(void)
    251251{
    252252    static bool threaded = false;       // Are we running threaded?

  • trunk/ppStack/src/ppStackThread.h

    r30620 r31158  
    4343
    4444// Initialise the threads
    45 void ppStackThreadInit(void);
     45void ppStackSetThreads(void);
    4646
    4747

  • trunk/ppStack/src/ppStackUpdateHeader.c

    r30620 r31158  
    7070    psMetadataAddStr(hdu->header, PS_LIST_TAIL, "TESS_ID", PS_META_REPLACE, "type of stack", tessID);
    7171
     72    psMetadataAddF32(hdu->header, PS_LIST_TAIL, "AIRM_SLP", PS_META_REPLACE, "airmass slope", options->airmassSlope);
     73   
     74
     75    // write the following fields to the output headers:
     76    // INP_NNNN : input file names
     77    // SCL_NNNN : scale factor for each input file
     78    // ZPT_NNNN : zero point of input file
     79    // EXT_NNNN : exptime of input file
     80    // AIR_NNNN : airmass of input file
     81
     82    // use separate loops so they are blocked together in the headers (more readable)
     83
     84    // save the input filenames
     85    for (int i = 0; i < options->num; i++) {
     86        char field[64];
     87        snprintf (field, 64, "INP_%04d", i);
     88
     89        psString basename = psStringFileBasename(options->origImages->data[i]);
     90        psMetadataAddStr(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image name", basename);
     91        psFree(basename);
     92    }   
     93
     94    // save the input normalizations
     95    for (int i = 0; i < options->num; i++) {
     96        char field[64];
     97
     98        float value = options->norm ? pow(10.0, -0.4*options->norm->data.F32[i]) : NAN;
     99        snprintf (field, 64, "SCL_%04d", i);
     100        psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image scale factor", value);
     101    }   
     102    // save the input exptimes
     103    for (int i = 0; i < options->num; i++) {
     104        char field[64];
     105
     106        float value = options->zpInput ? options->zpInput->data.F32[i] : NAN;
     107        snprintf (field, 64, "ZPT_%04d", i);
     108        psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image zero point", value);
     109    }   
     110    // save the input normalizations
     111    for (int i = 0; i < options->num; i++) {
     112        char field[64];
     113
     114        float value = options->expTimeInput ? options->expTimeInput->data.F32[i] : NAN;
     115        snprintf (field, 64, "EXP_%04d", i);
     116        psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image exptime", value);
     117    }   
     118    // save the input normalizations
     119    for (int i = 0; i < options->num; i++) {
     120        char field[64];
     121
     122        float value = options->airmassInput ? options->airmassInput->data.F32[i] : NAN;
     123        snprintf (field, 64, "AIR_%04d", i);
     124        psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image airmass", value);
     125    }   
    72126    return true;
    73127}
Note: See TracChangeset for help on using the changeset viewer.