IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 23704 for branches/pap/ppSub


Ignore:
Timestamp:
Apr 3, 2009, 11:58:24 AM (17 years ago)
Author:
Paul Price
Message:

Committing changes in progress. Attempting to make ppSub spit out not just the regular subtracted image and sources, but the inverse with sources in the inverse (for the warp-warp difference). Does not compile yet.

Location:
branches/pap/ppSub
Files:
4 added
10 edited

Legend:

Unmodified
Added
Removed

  • branches/pap/ppSub/configure.ac

    r21038 r23704  
    2323PKG_CHECK_MODULES([PSPHOT], [psphot >= 0.9.0])
    2424
     25AC_PATH_PROG([ERRORCODES], [psParseErrorCodes], [missing])
     26if test "$ERRORCODES" = "missing" ; then
     27  AC_MSG_ERROR([psParseErrorCodes is required])
     28fi
     29
    2530dnl Set CFLAGS for build
    2631IPP_STDOPTS

  • branches/pap/ppSub/src/Makefile.am

    r23688 r23704  
    4242        ppSub.h
    4343
     44
     45### Error codes.
     46BUILT_SOURCES = ppSubErrorCodes.h ppSubErrorCodes.c
     47CLEANFILES = ppSubErrorCodes.h ppSubErrorCodes.c
     48
     49ppSubErrorCodes.h : ppSubErrorCodes.dat ppSubErrorCodes.h.in
     50        $(ERRORCODES) --data=ppSubErrorCodes.dat --outdir=. ppSubErrorCodes.h
     51
     52ppSubErrorCodes.c : ppSubErrorCodes.dat ppSubErrorCodes.c.in ppSubErrorCodes.h
     53        $(ERRORCODES) --data=ppSubErrorCodes.dat --outdir=. ppSubErrorCodes.c
     54
     55
    4456clean-local:
    4557        -rm -f TAGS
    … …  
    4860tags:
    4961        etags `find . -name \*.[ch] -print`
    50 

  • branches/pap/ppSub/src/ppSub.c

    r23242 r23704  
    4949    }
    5050
    51     if (!ppSubArgumentsSetup(argc, argv, config)) {
     51    ppSubData *data = ppSubDataAlloc(); // Processing data
     52
     53    if (!ppSubArgumentsSetup(argc, argv, config, data)) {
    5254        psErrorStackPrint(stderr, "Error reading arguments.\n");
    5355        exitValue = PS_EXIT_CONFIG_ERROR;
    … …  
    5557    }
    5658
    57     if (!ppSubCamera(config)) {
     59    if (!ppSubCamera(config, data)) {
    5860        psErrorStackPrint(stderr, "Error setting up camera.\n");
    5961        exitValue = PS_EXIT_CONFIG_ERROR;
    … …  
    6163    }
    6264
    63     if (!ppSubArgumentsParse(config)) {
    64         psErrorStackPrint(stderr, "Error reading arguments.\n");
    65         exitValue = PS_EXIT_CONFIG_ERROR;
    66         goto die;
    67     }
    68 
    69     if (!ppSubLoop(config)) {
     65    if (!ppSubLoop(config, data)) {
    7066        psErrorStackPrint(stderr, "Error performing subtraction.\n");
    7167        exitValue = PS_EXIT_SYS_ERROR;

  • branches/pap/ppSub/src/ppSub.h

    r23688 r23704  
    2424// Output files, for activation/deactivation
    2525typedef enum {
    26     PPSUB_FILES_IMAGE = 0x01,           // Image files
    27     PPSUB_FILES_PHOT  = 0x02,           // Photometry files
    28     PPSUB_FILES_ALL   = 0xFF,           // All files
     26    PPSUB_FILES_INPUT    = 0x01,         // Input files
     27    PPSUB_FILES_CONV     = 0x02,         // Convolved files (output)
     28    PPSUB_FILES_SUB      = 0x04,         // Subtracted files (output)
     29    PPSUB_FILES_INV      = 0x08,         // Inverse subtracted files (output)
     30    PPSUB_FILES_PHOT_SUB = 0x10,         // Photometry files (output)
     31    PPSUB_FILES_PHOT_INV = 0x20,        // Photometry files (
     32    PPSUB_FILES_ALL      = 0xFF,         // All files
    2933} ppSubFiles;
    3034
    3135/// Data for processing
    3236typedef struct {
    33     psErrorCode quality;                /// Quality code; 0 for no problem
    34     psMetadata *stats;                  /// Statistics
     37    psErrorCode quality;                // Quality code; 0 for no problem
     38    bool photometry;                    // Perform photometry?
     39    bool inverse;                       // Output inverse subtraction as well?
     40    const char *statsFile;              // Statistics file
     41    psMetadata *stats;                  // Statistics
    3542} ppSubData;
    3643

  • branches/pap/ppSub/src/ppSubArguments.c

    r23402 r23704  
    189189}
    190190
    191 bool ppSubArgumentsSetup(int argc, char *argv[], pmConfig *config)
    192 {
    193     //    int argnum;                         // Argument Number
     191bool ppSubArguments(int argc, char *argv[], pmConfig *config, ppSubData *data)
     192{
    194193    assert(config);
    195 
    196     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    197     if (!recipe) {
    198         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSUB_RECIPE);
    199         return false;
    200     }
    201 
    202194
    203195    psMetadata *arguments = config->arguments; // Command-line arguments
    … …  
    212204    psMetadataAddStr(arguments, PS_LIST_TAIL, "-kernel", 0, "Precalculated kernel to apply", NULL);
    213205    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL);
    214     psMetadataAddF32(arguments, PS_LIST_TAIL, "-region", 0, "Size of iso-kernel region", NAN);
    215     psMetadataAddS32(arguments, PS_LIST_TAIL, "-size", 0, "Kernel half-size (pixels)", 0);
    216     psMetadataAddS32(arguments, PS_LIST_TAIL, "-order", 0, "Spatial polynomial order", -1);
    217     psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0,
    218                      "Kernel type (ISIS|POIS|SPAM|FRIES|GUNK|RINGS)", NULL);
    219     psMetadataAddF32(arguments, PS_LIST_TAIL, "-penalty", 0, "Penalty for wideness", NAN);
    220     psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-widths", 0,
    221                      "ISIS Gaussian FWHMs (comma-separated)", NULL);
    222     psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-orders", 0,
    223                      "ISIS polynomial orders (comma-separated)", NULL);
    224     psMetadataAddS32(arguments, PS_LIST_TAIL, "-rings-order", 0, "RINGS polynomial order", -1);
    225     psMetadataAddS32(arguments, PS_LIST_TAIL, "-inner", 0, "SPAM and FRIES inner radius", -1);
    226     psMetadataAddS32(arguments, PS_LIST_TAIL, "-spam-binning", 0, "SPAM kernel binning", 0);
    227     psMetadataAddF32(arguments, PS_LIST_TAIL, "-spacing", 0, "Typical stamp spacing (pixels)", NAN);
    228     psMetadataAddS32(arguments, PS_LIST_TAIL, "-footprint", 0, "Stamp footprint half-size (pixels)", -1);
    229     psMetadataAddF32(arguments, PS_LIST_TAIL, "-source-radius", 0, "Source matching radius (pixels)", NAN);
    230     psMetadataAddS32(arguments, PS_LIST_TAIL, "-stride", 0, "Size of convolution patches (pixels)", -1);
    231     psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold", 0, "Minimum threshold for stamps (ADU)", NAN);
    232     psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", -1);
    233     psMetadataAddF32(arguments, PS_LIST_TAIL, "-rej", 0, "Rejection thresold (sigma)", NAN);
    234     psMetadataAddF32(arguments, PS_LIST_TAIL, "-sys", 0, "Relative systematic error in kernel", NAN);
    235     psMetadataAddImageMask(arguments,  PS_LIST_TAIL, "-mask-bad", 0, "Mask value for bad pixels", 0);
    236     psMetadataAddImageMask(arguments,  PS_LIST_TAIL, "-mask-poor", 0, "Mask value for poor pixels", 0);
    237     psMetadataAddF32(arguments,  PS_LIST_TAIL, "-poor-frac", 0, "Fraction of variance for poor pixels", NAN);
    238     psMetadataAddF32(arguments, PS_LIST_TAIL, "-badfrac", 0, "Maximum fraction of bad pixels to accept", 1.0);
    239     psMetadataAddBool(arguments,  PS_LIST_TAIL, "-reverse", 0, "Reverse sense of subtraction?", false);
    240     psMetadataAddBool(arguments,  PS_LIST_TAIL, "-generate-mask", 0, "Generate mask if not supplied?", false);
    241     psMetadataAddStr(arguments,  PS_LIST_TAIL, "-stamps", 0,
    242                      "Stamps filename; file has x,y on each line", NULL);
    243     psMetadataAddBool(arguments, PS_LIST_TAIL, "-opt", 0,
    244                       "Derive optimum parameters for ISIS kernels?", false);
    245     psMetadataAddF32(arguments, PS_LIST_TAIL, "-opt-min", 0, "Minimum value for optimum kernel search", NAN);
    246     psMetadataAddF32(arguments, PS_LIST_TAIL, "-opt-max", 0, "Minimum value for optimum kernel search", NAN);
    247     psMetadataAddF32(arguments, PS_LIST_TAIL, "-opt-step", 0, "Step value for optimum kernel search", NAN);
    248     psMetadataAddF32(arguments, PS_LIST_TAIL, "-opt-tol", 0, "Tolerance for optimum kernel search", NAN);
    249     psMetadataAddS32(arguments, PS_LIST_TAIL, "-opt-order", 0, "Maximum order for optimum kernel search", -1);
    250     psMetadataAddBool(arguments, PS_LIST_TAIL, "-dual", 0, "Dual convolution", false);
    251     psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Perform photometry?", false);
     206    psMetadataAddStr(arguments,  PS_LIST_TAIL, "-stamps", 0, "Stamps filename; x,y on each line", NULL);
    252207    psMetadataAddS32(arguments, PS_LIST_TAIL, "-threads", 0, "Number of threads", 0);
    253     psMetadataAddS32(arguments, PS_LIST_TAIL, "-bin1", 0, "Binning factor for first level", 0);
    254     psMetadataAddS32(arguments, PS_LIST_TAIL, "-bin2", 0, "Binning factor for second level", 0);
    255208    psMetadataAddStr(arguments, PS_LIST_TAIL, "-dumpconfig", 0, "file to dump configuration to", NULL);
    256209    psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "Show diagnostic plots", NULL);
    … …  
    302255    }
    303256
    304     if (psMetadataLookupBool(NULL, arguments, "-photometry")) {
    305         psMetadataAddBool(recipe, PS_LIST_TAIL, "PHOTOMETRY", PS_META_REPLACE, "Perform photometry?", true);
    306     }
    307 
    308     return true;
    309 }
    310 
    311 bool ppSubArgumentsParse(pmConfig *config)
    312 {
    313     assert(config);
    314 
    315     psMetadata *arguments = config->arguments; // Command-line arguments
    316 
    317     valueArgStr(arguments, "-stats",  "STATS",  arguments);
    318     valueArgStr(arguments, "-stamps", "STAMPS", arguments);
    319 
    320     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    321     if (!recipe) {
    322         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSUB_RECIPE);
    323         goto ERROR;
    324     }
    325 
    326     VALUE_ARG_RECIPE_FLOAT("-region",        "REGION.SIZE",     F32);
    327     VALUE_ARG_RECIPE_INT("-size",            "KERNEL.SIZE",     S32, 0);
    328     VALUE_ARG_RECIPE_INT("-order",           "SPATIAL.ORDER",   S32, -1);
    329     VALUE_ARG_RECIPE_FLOAT("-spacing",       "STAMP.SPACING",   F32);
    330     VALUE_ARG_RECIPE_INT("-rings-order",     "RINGS.ORDER",     S32, -1);
    331     VALUE_ARG_RECIPE_INT("-inner",           "INNER",           S32, -1);
    332     VALUE_ARG_RECIPE_INT("-spam-binning",    "SPAM.BINNING",    S32, 0);
    333     VALUE_ARG_RECIPE_INT("-footprint",       "STAMP.FOOTPRINT", S32, -1);
    334     VALUE_ARG_RECIPE_FLOAT("-source-radius", "SOURCE.RADIUS",   F32);
    335     VALUE_ARG_RECIPE_INT("-stride",          "STRIDE",          S32, -1);
    336     VALUE_ARG_RECIPE_FLOAT("-threshold",     "STAMP.THRESHOLD", F32);
    337     VALUE_ARG_RECIPE_INT("-iter",            "ITER",            S32, -1);
    338     VALUE_ARG_RECIPE_FLOAT("-rej",           "REJ",             F32);
    339     VALUE_ARG_RECIPE_FLOAT("-sys",           "SYS",             F32);
    340     VALUE_ARG_RECIPE_FLOAT("-badfrac",       "BADFRAC",         F32);
    341     VALUE_ARG_RECIPE_FLOAT("-penalty",       "PENALTY",         F32);
    342     VALUE_ARG_RECIPE_FLOAT("-poor-frac",     "POOR.FRACTION",   F32);
    343     VALUE_ARG_RECIPE_INT("-bin1",            "BIN1",            S32, 0);
    344     VALUE_ARG_RECIPE_INT("-bin2",            "BIN2",            S32, 0);
    345 
    346     valueArgRecipeStr(arguments, recipe, "-mask-in",   "MASK.IN",  recipe);
    347     valueArgRecipeStr(arguments, recipe, "-mask-bad",  "MASK.BAD",  recipe);
    348     valueArgRecipeStr(arguments, recipe, "-mask-poor", "MASK.POOR", recipe);
    349 
    350     vectorArgRecipe(arguments, "-isis-widths", recipe, "ISIS.WIDTHS", recipe, PS_TYPE_F32);
    351     vectorArgRecipe(arguments, "-isis-orders", recipe, "ISIS.ORDERS", recipe, PS_TYPE_S32);
    352 
    353     psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
    354     psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders
    355     if (widths->n != orders->n) {
    356         psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Size of vectors for ISIS widths and orders do not match.");
    357         goto ERROR;
    358     }
    359 
    360     if (psMetadataLookupBool(NULL, arguments, "-opt") || psMetadataLookupBool(NULL, recipe, "OPTIMUM")) {
    361         psMetadataAddBool(recipe, PS_LIST_TAIL, "OPTIMUM", PS_META_REPLACE,
    362                           "Derive optimum parameters for ISIS kernels?", true);
    363         VALUE_ARG_RECIPE_FLOAT("-opt-min", "OPTIMUM.MIN",   F32);
    364         VALUE_ARG_RECIPE_FLOAT("-opt-max", "OPTIMUM.MAX",   F32);
    365         VALUE_ARG_RECIPE_FLOAT("-opt-step","OPTIMUM.STEP",  F32);
    366         VALUE_ARG_RECIPE_FLOAT("-opt-tol", "OPTIMUM.TOL",   F32);
    367         VALUE_ARG_RECIPE_INT("-opt-order", "OPTIMUM.ORDER", S32, -1);
    368     }
    369 
    370     psMetadataAddBool(arguments, PS_LIST_TAIL, "REVERSE", 0, "Reverse sense of subtraction",
    371                       psMetadataLookupBool(NULL, arguments, "-reverse"));
    372     psMetadataAddBool(recipe, PS_LIST_TAIL, "MASK.GENERATE", PS_META_REPLACE, "Generate mask if not supplied",
    373                       psMetadataLookupBool(NULL, arguments, "-generate-mask"));
    374     psMetadataAddBool(recipe, PS_LIST_TAIL, "DUAL", PS_META_REPLACE, "Dual convolution?",
    375                       psMetadataLookupBool(NULL, arguments, "-dual"));
    376 
    377     // Need to update this because it could have been overwritten by the camera's own recipe
    378     if (psMetadataLookupBool(NULL, arguments, "-photometry")) {
    379         psMetadataAddBool(recipe, PS_LIST_TAIL, "PHOTOMETRY", PS_META_REPLACE, "Perform photometry?", true);
     257    data->stamps = psMetadataLookupStr(NULL, arguments, "-stamps");
     258
     259    const char *statsName = psMetadataLookupStr(NULL, arguments, "-stats"); // Filename for statistics
     260    if (statsName && strlen(statsName) > 0) {
     261        psString resolved = pmConfigConvertFilename(statsName, config, true, true); // Resolved filename
     262        data->statsFile = fopen(resolved, "w");
     263        if (!data->statsFile) {
     264            psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
     265            psFree(resolved);
     266            return false;
     267        }
     268        psFree(resolved);
    380269    }
    381270
    … …  
    384273    }
    385274
    386     // Translate the kernel type
    387     psString type = psMetadataLookupStr(NULL, arguments, "-type"); // Name of kernel type
    388     if (!type || strlen(type) == 0) {
    389         type = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE");
    390         if (!type || strlen(type) == 0) {
    391             psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find KERNEL.TYPE specified.");
    392             goto ERROR;
    393         }
    394     }
    395     psMetadataAddStr(recipe, PS_LIST_TAIL, "KERNEL.TYPE", PS_META_REPLACE, "Type of kernel", type);
    396 
    397     psTrace("ppSub", 1, "Done reading command-line arguments\n");
    398 
    399     // XXX move this to ppSubArguments
    400     int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
     275    int threads = psMetadataLookupS32(NULL, arguments, "-threads"); // Number of threads
    401276    if (threads > 0) {
    402277        if (!psThreadPoolInit(threads)) {
    … …  
    407282
    408283    return true;
    409 
    410 ERROR:
    411     return false;
    412 }
     284}
     285
     286
     287
     288    psTrace("ppSub", 1, "Done reading command-line arguments\n");
     289
     290    return true;
     291}

  • branches/pap/ppSub/src/ppSubCamera.c

    r23688 r23704  
    161161
    162162
     163    // Now that the camera has been determined, we can read the recipe
     164    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
     165    if (!recipe) {
     166        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSUB_RECIPE);
     167        return false;
     168    }
     169    data->invert = psMetadataLookupBool(NULL, recipe, "INVERSE");
     170    data->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY");
     171
     172
    163173    // Output image
    164174    pmFPAfile *output = defineOutputFile(config, input, true, "PPSUB.OUTPUT", PM_FPA_FILE_IMAGE);
    … …  
    170180    output->save = true;
    171181    outMask->save = true;
    172     pmFPAfile *outVar = NULL;
    173182    if (inVar && refVar) {
    174         outVar = defineOutputFile(config, output, false, "PPSUB.OUTPUT.VARIANCE", PM_FPA_FILE_VARIANCE);
     183        pmFPAfile *outVar = defineOutputFile(config, output, false, "PPSUB.OUTPUT.VARIANCE",
     184                                             PM_FPA_FILE_VARIANCE);
    175185        if (!outVar) {
    176186            psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
    … …  
    180190    }
    181191
     192    if (data->inverse) {
     193        // Inverse output image
     194        pmFPAfile *outinv = defineOutputFile(config, input, true, "PPSUB.OUTINV", PM_FPA_FILE_IMAGE);
     195        pmFPAfile *outinvMask = defineOutputFile(config, outinv, false, "PPSUB.OUTINV.MASK",
     196                                                 PM_FPA_FILE_MASK);
     197        if (!outinv || !outinvMask) {
     198            psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
     199            return false;
     200        }
     201        outinv->save = true;
     202        outinvMask->save = true;
     203        if (inVar && refVar) {
     204            pmFPAfile *outinvVar = defineOutputFile(config, outinv, false, "PPSUB.OUTINV.VARIANCE",
     205                                                    PM_FPA_FILE_VARIANCE);
     206            if (!outinvVar) {
     207                psError(PS_ERR_UNKNOWN, false, "Unable to define output files");
     208                return false;
     209            }
     210            outinvVar->save = true;
     211        }
     212    }
    182213
    183214    // Convolved input image
    … …  
    252283
    253284    // psPhot input
    254     if (psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
     285    if (data->photometry) {
    255286        psphotModelClassInit();        // load implementation-specific models
    256287

  • branches/pap/ppSub/src/ppSubData.c

    r23688 r23704  
    2020
    2121
    22 static void subOptionsFree(ppSubData *options)
     22static void subDataFree(ppSubData *data)
    2323{
    24     psFree(options->stats);
     24    if (data->statsFile) {
     25        psString stats = psMetadataConfigFormat(data->stats); // Statistics to output
     26        if (!stats || strlen(stats) == 0) {
     27            psWarning("Unable to generate statistics file.");
     28        } else {
     29            fprintf(data->statsFile, "%s", stats);
     30        }
     31        psFree(stats);
     32        fclose(data->statsFile);
     33    }
     34    psFree(data->stats);
    2535    return;
    2636}
    … …  
    2838ppSubData *ppSubDataAlloc(void)
    2939{
    30     ppSubData *options = psAlloc(sizeof(ppSubData)); // Processing data, to return
    31     psMemSetDeallocator(options, (psFreeFunc)subOptionsFree);
     40    ppSubData *data = psAlloc(sizeof(ppSubData)); // Processing data, to return
     41    psMemSetDeallocator(data, (psFreeFunc)subDataFree);
    3242
    33     options->quality = 0;
    34     options->stats = psMetadataAlloc();
    35     psMetadataAddS32(options->stats, PS_LIST_TAIL, "QUALITY", 0, "Data quality", 0);
     43    data->quality = 0;
     44    data->photometry = false;
     45    data->inverse = false;
     46    data->statsFile = NULL;
     47    data->stats = psMetadataAlloc();
     48    psMetadataAddS32(data->stats, PS_LIST_TAIL, "QUALITY", 0, "Data quality", 0);
    3649
    3750    return options;

  • branches/pap/ppSub/src/ppSubDefineOutput.c

    r23403 r23704  
    2121#include "ppSub.h"
    2222
    23 bool ppSubDefineOutput(pmConfig *config, const pmFPAview *view)
     23bool ppSubDefineOutput(const char *name, pmConfig *config, const ppSubData *data, const pmFPAview *view)
    2424{
     25    psAssert(name, "Require name");
    2526    psAssert(config, "Require configuration");
    2627    psAssert(view, "Require view");
    2728
    28     pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
     29    pmCell *outCell = pmFPAfileThisCell(config->files, view, name); // Output cell
    2930    pmFPA *outFPA = outCell->parent->parent; // Output FPA
    3031    pmHDU *outHDU = outFPA->hdu; // Output HDU
    … …  
    3334    }
    3435
    35     // generate an output readout (first check if it's already there by virtue of kernels)
    36     pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel
    37     if (!outRO) {
    38         outRO = pmReadoutAlloc(outCell); // Output readout: subtraction
     36    pmReadout *outRO = NULL;            // Output readout
     37    if (outCell->readouts && outCell->readouts->n > 0 && outCell->readouts->data[0]) {
     38        outRO = psMemIncrRefCounter(outCell->readouts->data[0]);
     39    } else {
     40        outRO = pmReadoutAlloc(outCell);
    3941    }
    4042
    … …  
    5557    if (!kernels) {
    5658        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SUBTRACTION.KERNEL");
    57         psFree(inConv);
    58         psFree(refConv);
    5959        psFree(outRO);
    6060        return false;
    … …  
    6565    outRO->analysis = psMetadataCopy(outRO->analysis, analysis);
    6666
    67 #ifdef TESTING
    68     {
    69         psImage *kernelImage = psMetadataLookupPtr(&mdok, analysis,
    70                                                    "SUBTRACTION.KERNEL.IMAGE"); // Image of kernel
    71         psFits *fits = psFitsOpen("kernel.fits", "w");
    72         psFitsWriteImage(fits, NULL, kernelImage, 0, NULL);
    73         psFitsClose(fits);
    74     }
    75 #endif
     67    psFree(outRO);
    7668
    7769    return true;

  • branches/pap/ppSub/src/ppSubLoop.c

    r23688 r23704  
    2121#include "ppSub.h"
    2222
    23 bool ppSubLoop(pmConfig *config)
     23bool ppSubLoop(pmConfig *config, ppSubData *data)
    2424{
    2525    psAssert(config, "Require configuration.");
    … …  
    2828    pmConfigRecipesCull(config, "PPSUB,PPSTATS,PSPHOT,MASKS,JPEG");
    2929
    30     ppSubData *data = ppSubDataAlloc(); // Processing data
    3130
    32     bool mdok;                          // Status of MD lookup
    33     const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); // Filename for statistics
    34     FILE *statsFile = NULL;             // File stream for statistics
    35     if (statsName && strlen(statsName) > 0) {
    36         psString resolved = pmConfigConvertFilename(statsName, config, true, true); // Resolved filename
    37         statsFile = fopen(resolved, "w");
    38         if (!statsFile) {
    39             psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
    40             psFree(resolved);
     31    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT");
     32    pmFPAfile *reference = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF");
     33    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT");
     34    psAssert(input && reference && output, "Require files");
     35
     36    if (!ppSubFilesIterateDown(config, PPSUB_FILES_INPUT | PPSUB_FILES_CONV)) {
     37        psError(PPSUB_ERR_IO, false, "Unable to load files.");
     38        return false;
     39    }
     40
     41    psTimerStart("PPSUB_MATCH");
     42
     43    if (!ppSubSetMasks(config)) {
     44        psError(PS_ERR_UNKNOWN, false, "Unable to set masks.");
     45        return false;
     46    }
     47
     48    if (!ppSubMatchPSFs(config, data)) {
     49        psError(PS_ERR_UNKNOWN, false, "Unable to match PSFs.");
     50        return false;
     51    }
     52    if (data->quality) {
     53        // Can't do anything at all
     54        return true;
     55    }
     56
     57    // Close input files
     58    if (!ppSubFilesIterateUp(config, PPSUB_FILES_INPUT)) {
     59        psError(PPSUB_ERR_IO, false, "Unable to close input files.");
     60        return false;
     61    }
     62
     63    // Set up subtraction files
     64    if (!ppSubFilesIterateDown(config, PPSUB_FILES_SUB)) {
     65        psError(PPSUB_ERR_IO, false, "Unable to set up subtraction files.");
     66        return false;
     67    }
     68
     69    if (!ppSubDefineOutput("PPSUB.OUTPUT", config, data, view)) {
     70        psError(PS_ERR_UNKNOWN, false, "Unable to define output.");
     71        return false;
     72    }
     73
     74    if (data->inverse && !ppSubDefineOutput("PPSUB.OUTINV", config, data, view)) {
     75        psError(PS_ERR_UNKNOWN, false, "Unable to define inverse.");
     76        return false;
     77    }
     78
     79    if (!data->quality && !ppSubMakePSF("PPSUB.OUTPUT", "PPSUB.OUTINV", config, data, view)) {
     80        psError(PS_ERR_UNKNOWN, false, "Unable to generate PSF.");
     81        return false;
     82    }
     83
     84    // Close convolved files
     85    if (!ppSubFilesIterateUp(config, PPSUB_FILES_CONV)) {
     86        psError(PPSUB_ERR_IO, false, "Unable to close input files.");
     87        return false;
     88    }
     89
     90    if (!ppSubReadoutSubtract("PPSUB.OUTPUT", config, view)) {
     91        psError(PS_ERR_UNKNOWN, false, "Unable to subtract images.");
     92        return false;
     93    }
     94
     95    // Higher order background subtraction using psphot
     96    if (!ppSubBackground("PPSUB.OUTPUT", config, view)) {
     97        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     98        return false;
     99    }
     100
     101    if (!data->quality && data->!ppSubReadoutPhotometry("PPSUB.OUTPUT", config, data, view)) {
     102        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry.");
     103        return false;
     104    }
     105
     106    if (!ppSubReadoutUpdate("PPSUB.OUTPUT", config, data, view)) {
     107        psError(PS_ERR_UNKNOWN, false, "Unable to update.");
     108        return false;
     109    }
     110
     111    // Perform statistics on the cell
     112    if (statsFile) {
     113        pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT"); // Output file
     114        if (!output) {
     115            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PPSUB.OUTPUT.\n");
    41116            goto ERROR;
    42117        }
    43         psFree(resolved);
     118        psImageMaskType maskValue = pmConfigMaskGet("MASK.VALUE", config);
     119        ppStatsFPA(data->stats, output->fpa, view, maskValue, config);
    44120    }
    45121
    46     pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT");
    47     if (!input) {
    48         psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find input data!\n");
    49         goto ERROR;
    50     }
    51 
    52     pmFPAfile *reference = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF");
    53     if (!reference) {
    54         psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find reference data!\n");
    55         goto ERROR;
    56     }
    57 
    58     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT");
    59     if (!output) {
    60         psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!\n");
    61         goto ERROR;
    62     }
    63 
    64     pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
    65 
    66     // Iterate over the FPA hierarchy
    67     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    68         goto ERROR;
    69     }
    70 
    71     pmChip *inChip;                    // Input chip of interest
    72     while ((inChip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
    73         pmChip *refChip = pmFPAviewThisChip(view, reference->fpa); // Reference chip of interest
    74         if ((!inChip->file_exists && refChip->file_exists) ||
    75             (inChip->file_exists && !refChip->file_exists)) {
    76             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "FPA format discrepency between input and reference");
    77             psFree(view);
    78             goto ERROR;
     122    if (data->inverse) {
     123        // Set up inverse subtraction files
     124        if (!ppSubFilesIterateDown(config, PPSUB_FILES_INV)) {
     125            psError(PPSUB_ERR_IO, false, "Unable to set up inverse files.");
     126            return false;
    79127        }
    80128
    81         if (!inChip->file_exists) {
    82             continue;
     129        if (!ppSubReadoutInverse("PPSUB.OUTINV", "PPSUB.OUTPUT", config, view)) {
     130            psError(PS_ERR_UNKNOWN, false, "Unable to invert images.");
     131            return false;
    83132        }
    84133
    85         if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    86             goto ERROR;
     134        // Close subtraction files and open inverse photometry files
     135        if (!ppSubFilesIterateUp(config, PPSUB_FILES_SUB)) {
     136            psError(PPSUB_ERR_IO, false, "Unable to close subtraction files.");
     137            return false;
    87138        }
    88139
    89         pmCell *inCell;                // Cell of interest
    90         while ((inCell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
    91             pmCell *refCell = pmFPAviewThisCell(view, reference->fpa); // Reference cell of interest
    92             if ((!inCell->file_exists && refCell->file_exists) ||
    93                 (inCell->file_exists && !refCell->file_exists)) {
    94                 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    95                         "FPA format discrepency between input and reference");
    96                 psFree(view);
    97                 goto ERROR;
    98             }
    99             if (!inCell->file_exists) {
    100                 continue;
    101             }
    102             if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    103                 goto ERROR;
    104             }
    105 
    106             pmReadout *inRO;           // Readin of interest
    107             while ((inRO = pmFPAviewNextReadout(view, input->fpa, 1))) {
    108                 if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    109                     goto ERROR;
    110                 }
    111                 pmReadout *refRO = pmFPAviewThisReadout(view, reference->fpa);// Reference readout of interest
    112                 if (!refRO || (!inRO->data_exists && refRO->data_exists) ||
    113                     (inRO->data_exists && !refRO->data_exists)) {
    114                     psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    115                             "FPA format discrepency between input and reference");
    116                     psFree(view);
    117                     goto ERROR;
    118                 }
    119                 if (!inRO->data_exists) {
    120                     continue;
    121                 }
    122 
    123                 // Perform the analysis
    124                 if (!ppSubReadout(config, data, view)) {
    125                     psError(PS_ERR_UNKNOWN, false, "Unable to subtract images.\n");
    126                     goto ERROR;
    127                 }
    128 
    129                 if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    130                     goto ERROR;
    131                 }
    132             }
    133 
    134             // Perform statistics on the cell
    135             if (statsFile) {
    136                 pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT"); // Output file
    137                 if (!output) {
    138                     psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PPSUB.OUTPUT.\n");
    139                     goto ERROR;
    140                 }
    141                 psImageMaskType maskValue = pmConfigMaskGet("MASK.VALUE", config);
    142                 ppStatsFPA(data->stats, output->fpa, view, maskValue, config);
    143             }
    144 
    145             if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    146                 goto ERROR;
    147             }
     140        if (!data->quality && data->!ppSubReadoutPhotometry("PPSUB.OUTINV", config, data, view)) {
     141        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry.");
     142        return false;
    148143        }
    149 
    150         if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    151             goto ERROR;
     144    } else {
     145        // Close subtraction files
     146        if (!ppSubFilesIterateUp(config, PPSUB_FILES_SUB)) {
     147            psError(PPSUB_ERR_IO, false, "Unable to close subtraction files.");
     148            return false;
    152149        }
    153150    }
    154151
    155     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    156         goto ERROR;
     152
     153    if (!ppSubReadoutUpdate("PPSUB.OUTPUT", config, data, view)) {
     154        psError(PS_ERR_UNKNOWN, false, "Unable to update.");
     155        return false;
    157156    }
    158157
     158
    159159    psFree(view);
    160 
    161     // Write out summary statistics
    162     if (statsFile) {
    163         psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_SUB", 0, "Time for subtraction completion",
    164                          psTimerMark("ppSub"));
    165 
    166         const char *statsMDC = psMetadataConfigFormat(data->stats);
    167         if (!statsMDC || strlen(statsMDC) == 0) {
    168             psWarning("Unable to generate statistics MDC file.\n");
    169         } else {
    170             fprintf(statsFile, "%s", statsMDC);
    171         }
    172         psFree((void *)statsMDC);
    173         fclose(statsFile);
    174     }
    175160
    176161    psString dump_file = psMetadataLookupStr(&mdok, config->arguments, "-dumpconfig");
    … …  
    180165    }
    181166
    182     psFree(data);
    183167    return true;
    184168
    185169ERROR:
    186     psFree(data);
    187170    return false;
    188171}

  • branches/pap/ppSub/src/ppSubReadout.c

    r23688 r23704  
    2121#include "ppSub.h"
    2222
    23 bool ppSubReadout(pmConfig *config, ppSubData *data, const pmFPAview *view)
     23bool ppSubReadout(const char *name, bool reverse, pmConfig *config, ppSubData *data, const pmFPAview *view)
    2424{
    25     psTimerStart("PPSUB_MATCH");
     25    psAssert(name, "Require name");
     26    psAssert(config, "Require configuration");
     27    psAssert(data, "Require data");
     28    psAssert(view, "Require view");
    2629
    27     if (!ppSubSetMasks(config, view)) {
    28         psError(PS_ERR_UNKNOWN, false, "Unable to set masks.");
    29         return false;
    30     }
    31 
    32     if (!ppSubMatchPSFs(config, data, view)) {
    33         psError(PS_ERR_UNKNOWN, false, "Unable to match PSFs.");
    34         return false;
    35     } else if (data->quality) {
    36         // Can't do anything at all
    37         return true;
    38     }
    39 
    40     if (!ppSubDefineOutput(config, view)) {
     30    if (!ppSubDefineOutput(name, config, data, view)) {
    4131        psError(PS_ERR_UNKNOWN, false, "Unable to define output.");
    4232        return false;
    4333    }
    4434
    45     if (!data->quality && !ppSubMakePSF(config, data, view)) {
     35    if (!data->quality && !ppSubMakePSF(name, config, data, view)) {
    4636        psError(PS_ERR_UNKNOWN, false, "Unable to generate PSF.");
    4737        return false;
    4838    }
    4939
    50     if (!ppSubReadoutSubtract(config, view)) {
     40    if (!ppSubReadoutSubtract(name, reverse, config, view)) {
    5141        psError(PS_ERR_UNKNOWN, false, "Unable to subtract images.");
    5242        return false;
    … …  
    5444
    5545    // Higher order background subtraction using psphot
    56     if (!ppSubBackground(config, view)) {
     46    if (!ppSubBackground(name, config, view)) {
    5747        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    5848        return false;
    5949    }
    6050
    61     if (!data->quality && !ppSubReadoutPhotometry(config, data, view)) {
     51    if (!data->quality && data->!ppSubReadoutPhotometry(name, config, data, view)) {
    6252        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry.");
    6353        return false;
    6454    }
    6555
    66     if (!ppSubReadoutUpdate(config, data, view)) {
     56    if (!ppSubReadoutUpdate(name, config, data, view)) {
    6757        psError(PS_ERR_UNKNOWN, false, "Unable to update.");
    6858        return false;
    6959    }
    7060
     61
     62
     63
    7164    return true;
    7265}
Note: See TracChangeset for help on using the changeset viewer.