IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26898


Ignore:
Timestamp:
Feb 10, 2010, 7:41:23 PM (16 years ago)
Author:
eugene
Message:

updates from eam_branches/20091201

Location:
trunk/ppStack/src
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStack.h

    r26076 r26898  
    1313typedef enum {
    1414    PPSTACK_MASK_CAL    = 0x01,         // Photometric calibration failed
    15     PPSTACK_MASK_MATCH  = 0x02,         // PSF-matching failed
    16     PPSTACK_MASK_CHI2   = 0x04,         // Chi^2 too deviant
    17     PPSTACK_MASK_REJECT = 0x08,         // Rejection failed
    18     PPSTACK_MASK_BAD    = 0x10,         // Bad image (too many pixels rejected)
     15    PPSTACK_MASK_PSF    = 0x02,         // PSF measurement failed
     16    PPSTACK_MASK_MATCH  = 0x04,         // PSF-matching failed
     17    PPSTACK_MASK_CHI2   = 0x08,         // Chi^2 too deviant
     18    PPSTACK_MASK_REJECT = 0x10,         // Rejection failed
     19    PPSTACK_MASK_BAD    = 0x20,         // Bad image (too many pixels rejected)
    1920    PPSTACK_MASK_ALL    = 0xff          // All errors
    2021} ppStackMask;

  • trunk/ppStack/src/ppStackCamera.c

    r26076 r26898  
    373373        psMetadataLookupBool(&mdok, config->arguments, "-photometry"); // perform photometry
    374374    if (doPhotom) {
    375         // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by
    376         // psphotDefineFiles
     375        // This pmFPAfile, PSPHOT.INPUT, is just used as a carrier; output files (eg,
     376        // PSPHOT.RESID) are defined by psphotDefineFiles
    377377        pmFPAfile *psphotInput = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT");
    378378        if (!psphotInput) {
    … …  
    380380            return false;
    381381        }
     382        // specify the number of psphot input images
     383        psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", 1);
    382384
    383385        // Define associated psphot input/output files

  • trunk/ppStack/src/ppStackConvolve.c

    r26454 r26898  
    5959        pmFPAfileActivate(config->files, false, NULL);
    6060        ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i);
    61         if (options->convolve) {
    62             // XXX PPSTACK.CONV.KERNEL not defined unless convolve
    63             // pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL");
    64             pmFPAfileActivateSingle(config->files, true, "PPSTACK.CONV.KERNEL", i); // Activated file
    65         }
     61        if (options->convolve) {
     62            // XXX PPSTACK.CONV.KERNEL not defined unless convolve
     63            // pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL");
     64            pmFPAfileActivateSingle(config->files, true, "PPSTACK.CONV.KERNEL", i); // Activated file
     65        }
    6666
    6767        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
    … …  
    9494        options->origCovars->data[i] = psMemIncrRefCounter(readout->covariance);
    9595        if (!ppStackMatch(readout, options, i, config)) {
     96            // XXX many things can cause a failure of ppStackMatch -- should some be handled differently?
    9697            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
    9798            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH;

  • trunk/ppStack/src/ppStackFinish.c

    r26117 r26898  
    8383            fprintf(options->statsFile, "%s", statsMDC);
    8484        }
    85         psFree((void *)statsMDC);
     85        psFree(statsMDC);
    8686        fclose(options->statsFile); options->statsFile = NULL;
    8787        pmConfigRunFilenameAddWrite(config, "STATS", psMetadataLookupStr(NULL, config->arguments, "STATS"));

  • trunk/ppStack/src/ppStackMatch.c

    r26454 r26898  
    144144    psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
    145145
    146     psImage *binned = psphotBackgroundModel(ro, config); // Binned background model
     146    psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model
    147147    psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis,
    148148                                                  "PSPHOT.BACKGROUND.BINNING"); // Binning for model
    … …  
    274274                                                                PM_SUBTRACTION_ANALYSIS_KERNEL);
    275275
    276             pmSubtractionAnalysis(readout->analysis, NULL, kernels, region,
     276            pmSubtractionAnalysis(conv->analysis, NULL, kernels, region,
    277277                                  readout->image->numCols, readout->image->numRows);
    278278
    279279            psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
     280            bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
    280281            psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
     282            psImageCovarianceSetThreads(oldThreads);
    281283            psFree(readout->covariance);
    282284            readout->covariance = covar;
    … …  
    298300            float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
    299301            float kernelError = psMetadataLookupF32(NULL, ppsub, "KERNEL.ERR"); // Relative systematic error in kernel
     302            float normFrac = psMetadataLookupF32(NULL, ppsub, "NORM.FRAC"); // Fraction of window for normalisn windw
    300303            float sysError = psMetadataLookupF32(NULL, ppsub, "SYS.ERR"); // Relative systematic error in images
     304            float skyErr = psMetadataLookupF32(NULL, ppsub, "SKY.ERR"); // Additional error in sky
     305            float covarFrac = psMetadataLookupF32(NULL, ppsub, "COVAR.FRAC"); // Fraction for covariance calculation
     306
    301307            const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type
    302308            pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
    … …  
    315321            float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor"
    316322
     323            bool scale = psMetadataLookupBool(NULL, ppsub, "SCALE");        // Scale kernel parameters?
     324            float scaleRef = psMetadataLookupF32(NULL, ppsub, "SCALE.REF"); // Reference for scaling
     325            float scaleMin = psMetadataLookupF32(NULL, ppsub, "SCALE.MIN"); // Minimum for scaling
     326            float scaleMax = psMetadataLookupF32(NULL, ppsub, "SCALE.MAX"); // Maximum for scaling
     327            if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) {
     328                psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     329                        "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.",
     330                        scaleRef, scaleMin, scaleMax);
     331                return false;
     332            }
     333
     334
    317335            // These values are specified specifically for stacking
    318336            const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename
    … …  
    358376            fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK);
    359377
     378#if 1
    360379            // Add the background into the target image
    361380            psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
    362381            psBinaryOp(fake->image, fake->image, "+", bgImage);
    363382            psFree(bgImage);
     383#endif
    364384
    365385#ifdef TESTING
    … …  
    395415            if (kernel) {
    396416                if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis,
    397                                                stride, kernelError, maskVal, maskBad, maskPoor,
     417                                               stride, kernelError, covarFrac, maskVal, maskBad, maskPoor,
    398418                                               poorFrac, badFrac)) {
    399419                    psError(PS_ERR_UNKNOWN, false, "Unable to convolve images.");
    … …  
    408428                }
    409429            } else {
     430                // Scale the input parameters
     431                psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
     432                if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy,
     433                                                       options->inputSeeing->data.F32[index],
     434                                                       options->targetSeeing, scaleRef, scaleMin, scaleMax)) {
     435                    psError(PS_ERR_UNKNOWN, false, "Unable to scale kernel parameters");
     436                    psFree(fake);
     437                    psFree(optWidths);
     438                    psFree(stampSources);
     439                    psFree(conv);
     440                    psFree(widthsCopy);
     441                    if (threads > 0) {
     442                        pmSubtractionThreadsFinalize(readout, fake);
     443                    }
     444                    return false;
     445                }
     446
    410447                if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
    411                                         threshold, stampSources, stampsName, type, size, order, widths,
     448                                        threshold, stampSources, stampsName, type, size, order, widthsCopy,
    412449                                        orders, inner, ringsOrder, binning, penalty,
    413                                         optimum, optWidths, optOrder, optThresh, iter, rej, sysError,
    414                                         kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac,
    415                                         PM_SUBTRACTION_MODE_2)) {
     450                                        optimum, optWidths, optOrder, optThresh, iter, rej, normFrac,
     451                                        sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor,
     452                                        poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
    416453                    psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
    417454                    psFree(fake);
    … …  
    419456                    psFree(stampSources);
    420457                    psFree(conv);
     458                    psFree(widthsCopy);
    421459                    if (threads > 0) {
    422460                        pmSubtractionThreadsFinalize(readout, fake);
    … …  
    424462                    return false;
    425463                }
     464                psFree(widthsCopy);
    426465            }
    427466
    … …  
    495534            while ((item = psMetadataGetAndIncrement(iter))) {
    496535                assert(item->type == PS_DATA_UNKNOWN);
    497                 // Set the normalisation dimensions, since these will be otherwise unavailable when reading
    498                 // the images by scans.
    499536                pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
    500                 kernel->numCols = readout->image->numCols;
    501                 kernel->numRows = readout->image->numRows;
    502 
    503537                kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
    504538            }
    … …  
    526560        }
    527561
     562        // Kernel normalisation
     563        {
     564            double sum = 0.0;           // Sum of chi^2
     565            int num = 0;                // Number of measurements of chi^2
     566            psString regex = NULL;      // Regular expression
     567            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
     568            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
     569            psFree(regex);
     570            psMetadataItem *item = NULL;// Item from iteration
     571            while ((item = psMetadataGetAndIncrement(iter))) {
     572                assert(item->type == PS_TYPE_F32);
     573                float norm = item->data.F32; // Normalisation
     574                sum += norm;
     575                num++;
     576            }
     577            psFree(iter);
     578            float conv = sum/num;       // Mean normalisation from convolution
     579            float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
     580            float renorm =  stars / conv; // Renormalisation to apply
     581            psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n",
     582                     index, renorm, conv, stars);
     583            psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
     584            psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
     585        }
     586
    528587        // Reject image completely if the maximum deconvolution fraction exceeds the limit
    529588        float deconv = psMetadataLookupF32(NULL, conv->analysis,
    530589                                           PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
    531590        if (deconv > deconvLimit) {
    532             psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
    533                       deconv, deconvLimit);
     591            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n",
     592                      deconv, deconvLimit, index);
    534593            psFree(conv);
    535594            return NULL;

  • trunk/ppStack/src/ppStackOptions.c

    r26454 r26898  
    2020    psFree(options->convVariances);
    2121    psFree(options->psf);
     22    psFree(options->inputSeeing);
    2223    psFree(options->inputMask);
    2324    psFree(options->sourceLists);
    2425    psFree(options->norm);
     26    psFree(options->sources);
    2527    psFree(options->cells);
    2628    psFree(options->kernels);
    … …  
    4345    options->convolve = true;
    4446    options->matchZPs = true;
     47    options->photometry = false;
    4548    options->stats = NULL;
    4649    options->statsFile = NULL;
    … …  
    5457    options->num = 0;
    5558    options->psf = NULL;
     59    options->inputSeeing = NULL;
     60    options->targetSeeing = NAN;
    5661    options->inputMask = NULL;
    5762    options->sourceLists = NULL;
    5863    options->norm = NULL;
     64    options->sources = NULL;
    5965    options->cells = NULL;
    6066    options->kernels = NULL;

  • trunk/ppStack/src/ppStackOptions.h

    r26454 r26898  
    1010    bool convolve;                      // Convolve images?
    1111    bool matchZPs;                      // Adjust relative fluxes based on transparency analysis?
     12    bool photometry;                    // Perform photometry?
    1213    psMetadata *stats;                  // Statistics for output
    1314    FILE *statsFile;                    // File to which to write statistics
    … …  
    1819    // Prepare
    1920    pmPSF *psf;                         // Target PSF
     21    psVector *inputSeeing;              // Input seeing FWHMs
     22    float targetSeeing;                 // Target seeing FWHM
    2023    float sumExposure;                  // Sum of exposure times
    2124    psVector *inputMask;                // Mask for inputs
    2225    psArray *sourceLists;               // Individual lists of sources for matching
    2326    psVector *norm;                     // Normalisation for each image
     27    psArray *sources;                   // Matched sources
    2428    // Convolve
    2529    psArray *cells;                     // Cells for convolved images --- a handle for reading again

  • trunk/ppStack/src/ppStackPSF.c

    r25480 r26898  
    1414                  const psArray *psfs, const psVector *inputMask)
    1515{
     16    bool mdok = false;
     17
    1618#ifndef TESTING
    1719    // Get the recipe values
    … …  
    2426    int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF
    2527
     28    psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in
     29    if (!mdok || !maskValStr) {
     30        psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe");
     31        return false;
     32    }
     33    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask
     34
    2635    for (int i = 0; i < psfs->n; i++) {
    2736        if (inputMask->data.U8[i]) {
    … …  
    3241
    3342    // Solve for the target PSF
    34     pmPSF *psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel,
    35                                psfOrder, psfOrder);
     43    pmPSF *psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal);
    3644    if (!psf) {
    3745        psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");

  • trunk/ppStack/src/ppStackPhotometry.c

    r24216 r26898  
    1919    psAssert(recipe, "We've thrown an error on this before.");
    2020
    21     bool mdok;                          // Status of MD lookup
    22     if (!psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
     21    if (!options->photometry) {
    2322        // Nothing to do
    2423        return true;
    … …  
    6059                           "Bits to use for mark", markValue);
    6160
    62     pmCell *sourcesCell = pmFPAfileThisCell(config->files, photView, "PPSTACK.OUTPUT");
    63     psArray *inSources = psMetadataLookupPtr(NULL, sourcesCell->analysis, "PSPHOT.SOURCES");
     61    psArray *inSources = options->sources;
    6462    if (!inSources) {
    6563        psError(PS_ERR_UNKNOWN, false, "Unable to find input sources");
    … …  
    9290    if (options->stats) {
    9391        pmReadout *photRO = pmFPAviewThisReadout(photView, photFile->fpa); // Readout with the sources
    94         psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
    95         psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0,
    96                          "Number of sources detected", sources ? sources->n : 0);
    97         psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0,
    98                          "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
     92        pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // detections
     93        if (detections) {
     94            psAssert (detections->allSources, "missing sources?");
     95            psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", detections->allSources->n);
     96        } else {
     97            psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", 0);
     98        }
     99        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0, "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
    99100    }
    100101    psFree(photView);

  • trunk/ppStack/src/ppStackPrepare.c

    r26454 r26898  
    176176
    177177
    178         bool redoPhot = psMetadataLookupBool(NULL, recipe, "PHOT");
    179         psArray *sources = NULL;
    180 
    181         if (options->convolve || options->matchZPs || redoPhot) {
    182             pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
    183             sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
    184             if (!sources) {
    185                 psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
    186                 return NULL;
    187             }
    188             options->sourceLists->data[index] = psMemIncrRefCounter(sources);
    189         }
     178        bool redoPhot = psMetadataLookupBool(NULL, recipe, "PHOT");
     179
     180        pmDetections *detections = NULL;
     181        if (options->convolve || options->matchZPs || options->photometry || redoPhot) {
     182            pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
     183            detections = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.DETECTIONS"); // Sources
     184            if (!detections) {
     185                psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find source detections in readout.");
     186                return NULL;
     187            }
     188            psAssert (detections->allSources, "missing sources?");
     189
     190            options->sourceLists->data[index] = psMemIncrRefCounter(detections->allSources);
     191        }
    190192
    191193        // Re-do photometry if we don't trust the source lists
    … …  
    194196            pmFPAfileActivate(config->files, false, NULL);
    195197            ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index);
    196             if (options->convolve) {
    197                 pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL");
    198             }
     198            if (options->convolve) {
     199                pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL");
     200            }
    199201            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
    200202            pmFPAview *photView = ppStackFilesIterateDown(config);
    … …  
    206208            pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
    207209
    208             if (!ppStackInputPhotometer(ro, sources, config)) {
     210            if (!ppStackInputPhotometer(ro, detections->allSources, config)) {
    209211                psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
    210212                psFree(view);
    … …  
    225227    }
    226228    psFree(fileIter);
     229
     230    psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
     231    bool havePSFs = false;                                // Do we have any PSFs?
     232    options->inputSeeing = psVectorAlloc(num, PS_TYPE_F32);
     233    psVectorInit(options->inputSeeing, NAN);
     234    for (int i = 0; i < num; i++) {
     235        pmPSF *psf = psfs->data[i];     // PSF for image
     236        if (!psf) {
     237            continue;
     238        }
     239        havePSFs = true;
     240
     241        if (psf->trendNx > 0 && psf->trendNy > 0) {
     242            double sumFWHM = 0.0;           // FWHM for image
     243            int numFWHM = 0;                // Number of FWHM measurements
     244            for (int y = 0; y < psf->trendNy; y++) {
     245                float yPos = ((float)y + 0.5) / (float)psf->trendNy * numRows;
     246                for (int x = 0; x < psf->trendNx; x++) {
     247                    float xPos = ((float)x + 0.5) / (float)psf->trendNx * numCols;
     248                    float fwhm = pmPSFtoFWHM(psf, xPos, yPos); // FWHM for image
     249                    if (isfinite(fwhm)) {
     250                        sumFWHM += fwhm;
     251                        numFWHM++;
     252                    }
     253                }
     254            }
     255            if (numFWHM == 0) {
     256                options->inputSeeing->data.F32[i] = NAN;
     257                options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
     258                psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", i);
     259            } else {
     260                options->inputSeeing->data.F32[i] = sumFWHM / (float)numFWHM;
     261            }
     262        } else {
     263            options->inputSeeing->data.F32[i] = pmPSFtoFWHM(psf, 0.5 * numCols, 0.5 * numRows);
     264        }
     265
     266        psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]);
     267    }
     268    if (havePSFs) {
     269        psLogMsg("ppStack", PS_LOG_INFO, "%s", log);
     270    }
     271    psFree(log);
    227272
    228273    // Generate target PSF
    … …  
    237282        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN,
    238283                         "Target PSF for stack", options->psf);
     284        options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * numCols, 0.5 * numRows); // FWHM for target
     285        psLogMsg("ppStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);
    239286
    240287        pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.OUTPUT"); // Output chip

  • trunk/ppStack/src/ppStackReadout.c

    r26454 r26898  
    126126    int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
    127127
    128     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
     128    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
    129129    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    130130    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
    … …  
    219219    bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
    220220
    221     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
     221    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
    222222    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    223223    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits

  • trunk/ppStack/src/ppStackSetup.c

    r26454 r26898  
    2222    psAssert(recipe, "We've thrown an error on this before.");
    2323
    24     options->matchZPs = psMetadataLookupBool(NULL, recipe, "MATCH.ZERO.POINTS"); // Adjust zero points based on tranparency analysis?
     24    // XXX : switch to this name? options->matchZPs = psMetadataLookupBool(NULL, recipe, "MATCH.ZERO.POINTS"); // Adjust zero points based on tranparency analysis?
     25    options->matchZPs = psMetadataLookupBool(NULL, recipe, "ZP"); // Adjust zero points?
     26
     27    options->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); // Perform photometry?
    2528
    2629    options->convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images?

  • trunk/ppStack/src/ppStackSources.c

    r26454 r26898  
    6161    PS_ASSERT_PTR_NON_NULL(config, false);
    6262
    63     if (!options->matchZPs) {
    64         int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
    65         options->norm = psVectorAlloc(num, PS_TYPE_F32);
    66         psVectorInit (options->norm, 0.0);
    67 
    68         // XXX do I need to set this?
    69         // options->sumExposure = sumExpTime;
    70 
    71         return true;
     63    if (!options->matchZPs && !options->photometry) {
     64        int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
     65        options->norm = psVectorAlloc(num, PS_TYPE_F32);
     66        psVectorInit (options->norm, 0.0);
     67
     68        // XXX do I need to set this?
     69        // options->sumExposure = sumExpTime;
     70
     71        return true;
    7272    }
    7373
    … …  
    191191#endif
    192192
    193     psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
    194                                            iter2, starRej2, starSys2, starLimit,
    195                                            transIter, transRej, transThresh); // Transparencies for each image
    196     if (!trans) {
    197         psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
    198         return false;
    199     }
     193    if (options->matchZPs) {
     194        psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
     195                                               iter2, starRej2, starSys2, starLimit,
     196                                               transIter, transRej, transThresh); // Transparencies per image
     197        if (!trans) {
     198            psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
     199            return false;
     200        }
     201        for (int i = 0; i < trans->n; i++) {
     202            if (!isfinite(trans->data.F32[i])) {
     203                inputMask->data.U8[i] |= PPSTACK_MASK_CAL;
     204            }
     205        }
     206        // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
     207        // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
     208        // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
     209        // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
     210        // We don't need to know the magnitude zero point for the filter, since it cancels out
     211        if (options->matchZPs) {
     212            options->norm = psVectorAlloc(num, PS_TYPE_F32);
     213            for (int i = 0; i < num; i++) {
     214                if (!isfinite(trans->data.F32[i])) {
     215                    continue;
     216                }
     217                psArray *sources = sourceLists->data[i]; // Sources of interest
     218                float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];
     219                options->norm->data.F32[i] = magCorr;
     220                psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n",
     221                         i, magCorr);
     222
     223                for (int j = 0; j < sources->n; j++) {
     224                    pmSource *source = sources->data[j]; // Source of interest
     225                    if (!source) {
     226                        continue;
     227                    }
     228                    source->psfMag += magCorr;
     229                }
     230            }
     231        }
     232
     233#ifdef TESTING
     234        dumpMatches("source_mags.dat", num, matches, zp, trans);
     235#endif
     236        psFree(trans);
     237
     238#ifdef TESTING
     239        // Double check: all transparencies should be zero
     240        {
     241            psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
     242            if (!matches) {
     243                psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
     244                psFree(zp);
     245                return false;
     246            }
     247            psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
     248                                                   transThresh, starRej, starSys);
     249            for (int i = 0; i < num; i++) {
     250                fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
     251            }
     252            psFree(trans);
     253        }
     254#endif
     255    }
     256
    200257
    201258#if 0
    … …  
    216273#endif
    217274
    218 #ifdef TESTING
    219     dumpMatches("source_mags.dat", num, matches, zp, trans);
    220 #endif
    221 
    222     for (int i = 0; i < trans->n; i++) {
    223         if (!isfinite(trans->data.F32[i])) {
    224             inputMask->data.U8[i] |= PPSTACK_MASK_CAL;
    225         }
    226     }
    227 
    228     // Save best matches SOMEWHERE for future photometry
    229     // XXX this is a really poor output location; clean up the pmFPAfiles used in ppStack
    230     pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT");
    231     psArray *sourcesBest = psArrayAllocEmpty(matches->n);
    232 
    233     // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible
    234     int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required
    235     for (int i = 0; i < matches->n; i++) {
    236         pmSourceMatch *match = matches->data[i]; // Match of interest
    237         if (match->num < minMatches) {
    238             continue;
    239         }
    240 
    241         // We need to grab a single instance of this source: just take the first available
    242         int image = match->image->data.S32[0]; // Index of image
    243         int index = match->index->data.S32[0]; // Index of source within image
    244         psArray *sources = sourceLists->data[image]; // Sources for image
    245         pmSource *source = sources->data[index]; // Source of interest
    246 
    247         psArrayAdd(sourcesBest, sourcesBest->n, source);
    248     }
    249     psMetadataAdd(sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE,
    250                   "psphot sources", sourcesBest);
    251     psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
    252     psFree(sourcesBest);
     275    if (options->photometry) {
     276        // Save best matches for future photometry
     277        psArray *sourcesBest = options->sources = psArrayAllocEmpty(matches->n);
     278
     279        // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible
     280        int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required
     281        for (int i = 0; i < matches->n; i++) {
     282            pmSourceMatch *match = matches->data[i]; // Match of interest
     283            if (match->num < minMatches) {
     284                continue;
     285            }
     286
     287            // We need to grab a single instance of this source: just take the first available
     288            int image = match->image->data.S32[0]; // Index of image
     289            int index = match->index->data.S32[0]; // Index of source within image
     290            psArray *sources = sourceLists->data[image]; // Sources for image
     291            pmSource *source = sources->data[index]; // Source of interest
     292
     293            psArrayAdd(sourcesBest, sourcesBest->n, source);
     294        }
     295        psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
     296    }
    253297
    254298    psFree(matches);
    255 
    256     // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
    257     // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
    258     // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
    259     // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
    260     // We don't need to know the magnitude zero point for the filter, since it cancels out
    261     options->norm = psVectorAlloc(num, PS_TYPE_F32);
    262     for (int i = 0; i < num; i++) {
    263         if (!isfinite(trans->data.F32[i])) {
    264             continue;
    265         }
    266         psArray *sources = sourceLists->data[i]; // Sources of interest
    267         float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];
    268         options->norm->data.F32[i] = magCorr;
    269         psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr);
    270 
    271         for (int j = 0; j < sources->n; j++) {
    272             pmSource *source = sources->data[j]; // Source of interest
    273             if (!source) {
    274                 continue;
    275             }
    276             source->psfMag += magCorr;
    277         }
    278     }
    279     psFree(trans);
    280 
    281 #ifdef TESTING
    282     // Double check: all transparencies should be zero
    283     {
    284         psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
    285         if (!matches) {
    286             psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
    287             psFree(zp);
    288             return false;
    289         }
    290         psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
    291                                                transThresh, starRej, starSys);
    292         for (int i = 0; i < num; i++) {
    293             fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
    294         }
    295         psFree(trans);
    296     }
    297 #endif
    298299
    299300    return true;
Note: See TracChangeset for help on using the changeset viewer.