IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 13489


Ignore:
Timestamp:
May 23, 2007, 12:41:52 PM (19 years ago)
Author:
Paul Price
Message:

Bug fixes and small improvements to get it working properly. Seems to work OK now, but hasn't been tested rigourously.

Location:
trunk/ppStack/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStackArguments.c

    r13466 r13489  
    8080    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stat", 0, "Statistics file", NULL);
    8181    psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0);
    82     psMetadataAddF32(arguments, PS_LIST_TAIL, "-rej", 0, "Rejection thresold (sigma)", NAN);
     82    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-rej", 0, "Combination rejection thresold (sigma)", NAN);
     83    psMetadataAddF32(arguments, PS_LIST_TAIL, "-convolve-rej", 0, "Convolution rejection thresold (sigma)", NAN);
    8384    psMetadataAddF32(arguments, PS_LIST_TAIL, "-extent", 0, "Extent of convolution (sigma)", NAN);
    8485    psMetadataAddU8(arguments,  PS_LIST_TAIL, "-mask-bad", 0, "Mask value for bad pixels", 0);
    … …  
    8990    }
    9091
    91     psMetadataAddStr(config->arguments, PS_LIST_TAIL, "IMAGE.LIST", 0,
     92    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "IMAGES.LIST", 0,
    9293                     "Name of the input image list", argv[1]);
    9394    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "MASKS.LIST", 0,
    … …  
    104105    }
    105106
    106     VALUE_ARG_RECIPE_INT("-iter",       "ITER",       S32, 0);
    107     VALUE_ARG_RECIPE_FLOAT("-rej",      "REJ",        F32);
    108     VALUE_ARG_RECIPE_FLOAT("-extent",   "EXTENT",     F32);
    109     VALUE_ARG_RECIPE_INT("-mask-bad",   "MASK.BAD",   U8, 0);
    110     VALUE_ARG_RECIPE_INT("-mask-blank", "MASK.BLANK", U8, 0);
     107    VALUE_ARG_RECIPE_INT("-iter",           "ITER",         S32, 0);
     108    VALUE_ARG_RECIPE_FLOAT("-combine-rej",  "COMBINE.REJ",  F32);
     109    VALUE_ARG_RECIPE_FLOAT("-convolve-rej", "CONVOLVE.REJ", F32);
     110    VALUE_ARG_RECIPE_FLOAT("-extent",       "EXTENT",       F32);
     111    VALUE_ARG_RECIPE_INT("-mask-bad",       "MASK.BAD",     U8, 0);
     112    VALUE_ARG_RECIPE_INT("-mask-blank",     "MASK.BLANK",   U8, 0);
    111113
    112114    psTrace("ppStack", 1, "Done reading command-line arguments\n");

  • trunk/ppStack/src/ppStackCamera.c

    r13464 r13489  
    104104
    105105    // Output image
    106     pmFPAfile *output = pmFPAfileDefineOutput(config, NULL, "PPSTACK.OUTPUT");
     106    pmFPA *fpa = pmFPAConstruct(config->camera); // FPA to contain the output
     107    if (!fpa) {
     108        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration.");
     109        return false;
     110    }
     111    pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, "PPSTACK.OUTPUT");
     112    psFree(fpa);                        // Drop reference
    107113    if (!output) {
    108114        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT"));
    … …  
    113119        return false;
    114120    }
     121
     122    // Generate the correct structure
     123    pmFPALevel phuLevel = pmFPAPHULevel(output->format); // Level at which PHU goes
     124    pmFPAview *view = pmFPAviewAlloc(0);// View for current level
     125    if (phuLevel == PM_FPA_LEVEL_FPA) {
     126        if (!pmFPAAddSourceFromView(fpa, "Stack", view, output->format)) {
     127            psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
     128            psFree(fpa);
     129            psFree(view);
     130            return NULL;
     131        }
     132    } else {
     133        pmChip *chip;                       // Chip from FPA
     134        while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
     135            if (phuLevel == PM_FPA_LEVEL_CHIP) {
     136                if (!pmFPAAddSourceFromView(fpa, "Stack", view, output->format)) {
     137                    psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
     138                    psFree(fpa);
     139                    psFree(view);
     140                    return NULL;
     141                }
     142            } else {
     143                pmCell *cell;                   // Cell from chip
     144                while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
     145                    if (phuLevel == PM_FPA_LEVEL_CELL) {
     146                        if (!pmFPAAddSourceFromView(fpa, "Stack", view, output->format)) {
     147                            psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
     148                            psFree(fpa);
     149                            psFree(view);
     150                            return NULL;
     151                        }
     152                    }
     153                }
     154            }
     155        }
     156    }
     157    psFree(view);
    115158
    116159    // Output mask

  • trunk/ppStack/src/ppStackLoop.c

    r13464 r13489  
    3030    }
    3131
    32     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file
    33     if (!output) {
    34         psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!\n");
     32    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Token input file
     33    if (!input) {
     34        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find input data!\n");
    3535        return false;
    3636    }
    … …  
    4545
    4646    pmChip *chip;                       // Chip of interest
    47     while ((chip = pmFPAviewNextChip(view, output->fpa, 1)) != NULL) {
     47    while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
    4848        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    4949            return false;
    … …  
    5151
    5252        pmCell *cell;                // Cell of interest
    53         while ((cell = pmFPAviewNextCell(view, output->fpa, 1)) != NULL) {
     53        while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
    5454            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    5555                return false;
    … …  
    6767
    6868            pmReadout *readout;         // Readout of interest
    69             while ((readout = pmFPAviewNextReadout(view, output->fpa, 1))) {
     69            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1))) {
    7070                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    7171                    return false;
    … …  
    8585            // Perform statistics on the cell
    8686            if (stats) {
    87                 ppStats(stats, output->fpa, view, config);
     87                ppStats(stats, input->fpa, view, config);
    8888            }
    8989

  • trunk/ppStack/src/ppStackReadout.c

    r13464 r13489  
    1717bool ppStackReadout(pmConfig *config, const pmFPAview *view)
    1818{
    19     // XXX Somehow need to add an HDU in to the output so that we can actually write it out!
     19    // Get the recipe values
     20    int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations
     21    float combineRej = psMetadataLookupF32(NULL, config->arguments, "COMBINE.REJ"); // Combination threshold
     22    float convolveRej = psMetadataLookupF32(NULL, config->arguments, "CONVOLVE.REJ"); // Convolution threshold
     23    float extent = psMetadataLookupF32(NULL, config->arguments, "EXTENT"); // Extent of convolution
     24    psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
     25    psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
    2026
    2127    // Get the output target
    … …  
    3137                                                            "^PPSTACK.INPUT$"); // Iterator over input files
    3238    psMetadataItem *item;               // Item from iteration
     39    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
     40    psRandom *rng = psRandomAlloc(0, PS_RANDOM_TAUS); // Random number generator
    3341    while ((item = psMetadataGetAndIncrement(inputIter))) {
    3442        assert(item->type == PS_DATA_UNKNOWN);
    3543        pmFPAfile *inputFile = item->data.V; // An input file
    3644        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout
     45
     46        // Brain-dead background subtraction
     47        if (!psImageBackground(stats, ro->image, ro->mask, maskBad, rng)) {
     48            psError(PS_ERR_UNKNOWN, false, "Unable to get image background on image %ld", stack->n);
     49            psFree(stats);
     50            psFree(rng);
     51            psFree(inputIter);
     52            psFree(stack);
     53            psFree(outRO);
     54            return false;
     55        }
     56        (void)psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(stats->robustMedian, PS_TYPE_F32));
     57
     58        // Normalise each input by the exposure time
     59        float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time
     60        if (!isfinite(exposure)) {
     61            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     62                    "CELL.EXPOSURE is not set for input file %ld", stack->n);
     63            psFree(stats);
     64            psFree(rng);
     65            psFree(inputIter);
     66            psFree(outRO);
     67            psFree(stack);
     68            return false;
     69        }
     70        (void)psBinaryOp(ro->image, ro->image, "/", psScalarAlloc(exposure, PS_TYPE_F32));
    3771        pmStackData *data = pmStackDataAlloc(ro, SEEING, WEIGHT); // Data to stack
    3872        psArrayAdd(stack, ARRAY_BUFFER, data);
    3973    }
    4074    psFree(inputIter);
     75    psFree(stats);
     76    psFree(rng);
    4177
    42     // Get the recipe values
    43     int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations
    44     float rej = psMetadataLookupF32(NULL, config->arguments, "REJ"); // Rejection threshold
    45     float extent = psMetadataLookupF32(NULL, config->arguments, "EXTENT"); // Extent of convolution
    46     psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
    47     psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
    48 
    49     if (!pmStackCombine(outRO, stack, maskBad, maskBlank, iter, rej)) {
     78    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, iter, combineRej)) {
    5079        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
    5180        psFree(stack);
    … …  
    5483    }
    5584
    56     if (!pmStackReject(stack, maskBad, extent, rej)) {
     85#if 1
     86    {
     87        psFits *fits = psFitsOpen("combined.fits", "w");
     88        psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
     89        psFitsClose(fits);
     90    }
     91
     92    for (int i = 0; i < num; i++) {
     93        pmStackData *data = stack->data[i]; // Data for this image
     94        psImage *inspected = psPixelsToMask(NULL, data->pixels,
     95                                            psRegionSet(0, data->sky->image->numCols,
     96                                                        0, data->sky->image->numRows),
     97                                            maskBlank);
     98        psString name = NULL;           // Name of image
     99        psStringAppend(&name, "inspect%03d.fits", i);
     100        psFits *fits = psFitsOpen(name, "w");
     101        psFree(name);
     102        psFitsWriteImage(fits, NULL, inspected, 0, NULL);
     103        psFitsClose(fits);
     104    }
     105#endif
     106
     107    if (!pmStackReject(stack, maskBad, extent, convolveRej)) {
    57108        psError(PS_ERR_UNKNOWN, false, "Unable to reject input pixels.");
    58109        psFree(stack);
    … …  
    61112    }
    62113
    63     if (!pmStackCombine(outRO, stack, maskBad, maskBlank, 0, 0.0)) {
     114#if 1
     115    for (int i = 0; i < num; i++) {
     116        pmStackData *data = stack->data[i]; // Data for this image
     117        psImage *rejected = psPixelsToMask(NULL, data->pixels,
     118                                           psRegionSet(0, data->sky->image->numCols,
     119                                                       0, data->sky->image->numRows),
     120                                           maskBlank);
     121        psString name = NULL;           // Name of image
     122        psStringAppend(&name, "reject%03d.fits", i);
     123        psFits *fits = psFitsOpen(name, "w");
     124        psFree(name);
     125        psFitsWriteImage(fits, NULL, rejected, 0, NULL);
     126        psFitsClose(fits);
     127    }
     128#endif
     129
     130    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, 0, combineRej)) {
    64131        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
    65132        psFree(stack);
    … …  
    68135    }
    69136
     137    outRO->data_exists = true;
     138    outRO->parent->data_exists = true;
     139    outRO->parent->parent->data_exists = true;
     140
    70141    psFree(stack);
    71142    psFree(outRO);                      // Drop reference
Note: See TracChangeset for help on using the changeset viewer.