IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14107


Ignore:
Timestamp:
Jul 10, 2007, 2:15:40 PM (19 years ago)
Author:
Paul Price
Message:

Supporting new kernel types.

Location:
trunk/ppSub/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppSub/src/ppSubArguments.c

    r13750 r14107  
    1717{
    1818    fprintf(stderr, "\nPan-STARRS PSF-matched image subtraction\n\n");
    19     fprintf(stderr, "Usage: %s INPUT.fits REFERENCE.fits OUTPUT_ROOT \n"
    20             "\t[-psf REFERENCE.psf.fits|-psflist REFERENCE.list]\n",
    21             program);
     19    fprintf(stderr, "Usage: %s INPUT.fits REFERENCE.fits OUTPUT_ROOT\n", program);
    2220    fprintf(stderr, "\n");
    2321    psArgumentHelp(arguments);
    … …  
    5755    } \
    5856    psMetadataAdd##TYPE(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \
    59 }
    60 
    61 // Get a mask value from the command-line or recipe, and add it to the arguments
    62 #define VALUE_ARG_RECIPE_MASK(ARGNAME, RECIPENAME) { \
    63     bool mdok; \
    64     const char *name = psMetadataLookupStr(&mdok, arguments, ARGNAME); \
    65     if (!mdok || !name || strlen(name) == 0) { \
    66         name = psMetadataLookupStr(NULL, recipe, RECIPENAME); \
    67         if (!name) { \
    68             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \
    69                 RECIPENAME, PPSUB_RECIPE); \
    70             goto ERROR; \
    71         } \
    72     } \
    73     psMaskType value = pmConfigMask(name, config); \
    74     psMetadataAddU8(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \
    7557}
    7658
    … …  
    152134}
    153135
    154 // Add a single filename to the arguments as an array, so that it can be used with pmFPAfileBindFromArgs, etc
    155 static void fileList(const char *file, // The symbolic name for the file
    156                      const char *name, // The name of the file
    157                      const char *comment, // Description of the file
    158                      pmConfig *config // Configuration
    159     )
    160 {
    161     psArray *files = psArrayAlloc(1); // Array with file names
    162     files->data[0] = psStringCopy(name);
    163     psMetadataAddArray(config->arguments, PS_LIST_TAIL, file, 0, comment, files);
    164     psFree(files);
    165     return;
    166 }
    167 
    168 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
    169136
    170137bool ppSubArguments(int argc, char *argv[], pmConfig *config)
    171138{
    172139    assert(config);
    173 
    174     psMetadataAddBool(config->arguments, PS_LIST_TAIL, "PHOTOMETRY", 0, "Do photometry?",
    175                       (psArgumentGet(argc, argv, "-psf") || psArgumentGet(argc, argv, "-psflist")) ?
    176                       true : false);
    177     pmConfigFileSetsMD(config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist");
    178140
    179141    psMetadata *arguments = psMetadataAlloc(); // Command-line arguments
    … …  
    182144    psMetadataAddStr(arguments, PS_LIST_TAIL, "-refmask", 0, "Referencemask image", NULL);
    183145    psMetadataAddStr(arguments, PS_LIST_TAIL, "-refweight", 0, "Referenceweight image", NULL);
    184     psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL);
     146    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stat", 0, "Statistics file", NULL);
    185147
    186148    psMetadataAddS32(arguments, PS_LIST_TAIL, "-size", 0, "Kernel half-size (pixels)", 0);
    187149    psMetadataAddS32(arguments, PS_LIST_TAIL, "-order", 0, "Spatial polynomial order", 0);
    188     psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Kernel type (POIS or ISIS)", NULL);
     150    psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Kernel type (ISIS|POIS|SPAM|FRIES)", NULL);
    189151    psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-widths", 0, "ISIS Gaussian widths (comma-separated)", NULL);
    190152    psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-orders", 0, "ISIS polynomial orders (comma-separated)", NULL);
     153    psMetadataAddS32(arguments, PS_LIST_TAIL, "-inner", 0, "SPAM and FRIES inner radius", 0);
     154    psMetadataAddS32(arguments, PS_LIST_TAIL, "-spam-binning", 0, "SPAM kernel binning", 2);
    191155    psMetadataAddF32(arguments, PS_LIST_TAIL, "-spacing", 0, "Typical stamp spacing (pixels)", NAN);
    192156    psMetadataAddS32(arguments, PS_LIST_TAIL, "-footprint", 0, "Stamp footprint half-size (pixels)", 0);
    … …  
    202166    }
    203167
    204     fileList("INPUT", argv[1], "Name of the input image",     config);
    205     fileList("REF",   argv[2], "Name of the reference image", config);
     168    psArray *files = psArrayAlloc(1);   // Array with file names
     169    files->data[0] = psStringCopy(argv[1]);
     170    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "INPUT", 0, "Name of the input image", files);
     171    psFree(files);
     172    files = psArrayAlloc(1);
     173    files->data[0] = psStringCopy(argv[2]);
     174    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "REF", 0, "Name of the reference image", files);
     175    psFree(files);
    206176    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[3]);
    207177
    208     const char *inMask = psMetadataLookupStr(NULL, arguments, "-inmask"); // Name of input mask
    209     if (inMask && strlen(inMask) > 0) {
    210         fileList("INPUT.MASK", inMask, "Name of the input mask image", config);
    211     }
    212     const char *inWeight = psMetadataLookupStr(NULL, arguments, "-inweight"); // Name of input weight
    213     if (inWeight && strlen(inWeight) > 0) {
    214         fileList("INPUT.WEIGHT", inWeight, "Name of the input weight image", config);
    215     }
    216 
    217     const char *refMask = psMetadataLookupStr(NULL, arguments, "-refmask"); // Name of reference mask
    218     if (refMask && strlen(refMask) > 0) {
    219         fileList("REF.MASK", refMask, "Name of the reference mask image", config);
    220     }
    221     const char *refWeight = psMetadataLookupStr(NULL, arguments, "-refweight"); // Name of reference weight
    222     if (refWeight && strlen(refWeight) > 0) {
    223         fileList("REF.WEIGHT", refWeight, "Name of the reference weight image", config);
    224     }
    225 
    226     valueArgStr(config, arguments, "-stats",      "STATS",         config->arguments);
     178    valueArgStr(config, arguments, "-inmask",    "INPUT.MASK",    config->arguments);
     179    valueArgStr(config, arguments, "-inweight",  "INPUT.WEIGHT",  config->arguments);
     180    valueArgStr(config, arguments, "-refmask",   "REF.MASK",      config->arguments);
     181    valueArgStr(config, arguments, "-refweight", "REF.WEIGHT",    config->arguments);
     182    valueArgStr(config, arguments, "-stat",      "STATS",         config->arguments);
    227183
    228184    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    … …  
    232188    }
    233189
    234     VALUE_ARG_RECIPE_INT("-size",        "KERNEL.SIZE",     S32, 0);
    235     VALUE_ARG_RECIPE_INT("-order",       "SPATIAL.ORDER",   S32, 0);
    236     VALUE_ARG_RECIPE_FLOAT("-spacing",   "STAMP.SPACING",   F32);
    237     VALUE_ARG_RECIPE_INT("-footprint",   "STAMP.FOOTPRINT", S32, 0);
    238     VALUE_ARG_RECIPE_FLOAT("-threshold", "STAMP.THRESHOLD", F32);
    239     VALUE_ARG_RECIPE_INT("-iter",        "ITER",            S32, 0);
    240     VALUE_ARG_RECIPE_FLOAT("-rej",       "REJ",             F32);
    241     VALUE_ARG_RECIPE_MASK("-mask-bad",   "MASK.BAD");
    242     VALUE_ARG_RECIPE_MASK("-mask-blank", "MASK.BLANK");
     190    VALUE_ARG_RECIPE_INT("-size",         "KERNEL.SIZE",     S32, 0);
     191    VALUE_ARG_RECIPE_INT("-order",        "SPATIAL.ORDER",   S32, 0);
     192    VALUE_ARG_RECIPE_FLOAT("-spacing",    "STAMP.SPACING",   F32);
     193    VALUE_ARG_RECIPE_INT("-inner",        "INNER",           S32, 0);
     194    VALUE_ARG_RECIPE_INT("-spam-binning", "SPAM.BINNING",    S32, 0);
     195    VALUE_ARG_RECIPE_INT("-footprint",    "STAMP.FOOTPRINT", S32, 0);
     196    VALUE_ARG_RECIPE_FLOAT("-threshold",  "STAMP.THRESHOLD", F32);
     197    VALUE_ARG_RECIPE_INT("-iter",         "ITER",            S32, 0);
     198    VALUE_ARG_RECIPE_FLOAT("-rej",        "REJ",             F32);
     199    VALUE_ARG_RECIPE_INT("-mask-bad",     "MASK.BAD",        U8, 0);
     200    VALUE_ARG_RECIPE_INT("-mask-blank",   "MASK.BLANK",      U8, 0);
    243201
    244202    vectorArgRecipe(config, arguments, "-isis-widths", recipe, "ISIS.WIDTHS", config->arguments, PS_TYPE_F32);
    … …  
    269227    } else if (strcasecmp(type, "ISIS") == 0) {
    270228        kernelType = PM_SUBTRACTION_KERNEL_ISIS;
     229    } else if (strcasecmp(type, "SPAM") == 0) {
     230        kernelType = PM_SUBTRACTION_KERNEL_SPAM;
     231    } else if (strcasecmp(type, "FRIES") == 0) {
     232        kernelType = PM_SUBTRACTION_KERNEL_FRIES;
    271233    } else {
    272234        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", type);

  • trunk/ppSub/src/ppSubReadout.c

    r13738 r14107  
    66#include <pslib.h>
    77#include <psmodules.h>
    8 #include <psphot.h>
    98
    109#include "ppSub.h"
     10
     11#define MASK_BAD      0x01              // Mask value for bad pixel
     12#define MASK_STAMP    0x02              // Mask value for bad stamp (and bad stamp region)
    1113
    1214bool ppSubReadout(pmConfig *config, const pmFPAview *view)
    … …  
    1416    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
    1517    pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
     18#if 0
    1619    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
    1720    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
     21#endif
    1822
    1923    psImage *input = inRO->image;       // Input image
    … …  
    3438    psVector *widths = psMetadataLookupPtr(NULL, config->arguments, "ISIS.WIDTHS"); // ISIS Gaussian widths
    3539    psVector *orders = psMetadataLookupPtr(NULL, config->arguments, "ISIS.ORDERS"); // ISIS Polynomial orders
     40    int inner = psMetadataLookupS32(NULL, config->arguments, "INNER"); // Inner radius
     41    int binning = psMetadataLookupS32(NULL, config->arguments, "SPAM.BINNING"); // Binning for SPAM kernel
    3642    psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
    3743    psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
    3844
    39     int numCols = input->numCols, numRows = input->numRows; // Image dimensions
    40     if (!inRO->mask) {
    41         inRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
    42         psImageInit(inRO->mask, 0);
     45    if (!inRO->mask && !pmReadoutGenerateMask(inRO)) {
     46        psError(PS_ERR_UNKNOWN, false, "Unable to generate mask for input image");
     47        return false;
    4348    }
    44     if (!refRO->mask) {
    45         refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
    46         psImageInit(refRO->mask, 0);
     49    if (!inRO->weight && !pmReadoutGenerateWeight(inRO, true)) {
     50        psError(PS_ERR_UNKNOWN, false, "Unable to generate weight map for input image");
     51        return false;
     52    }
     53    if (!refRO->mask && !pmReadoutGenerateMask(refRO)) {
     54        psError(PS_ERR_UNKNOWN, false, "Unable to generate mask for reference image");
     55        return false;
     56    }
     57    if (!refRO->weight && !pmReadoutGenerateWeight(refRO, true)) {
     58        psError(PS_ERR_UNKNOWN, false, "Unable to generate weight map for reference image");
     59        return false;
    4760    }
    4861
    49     // Mask for subtraction
    50     psImage *subMask = pmSubtractionMask(inRO->mask, refRO->mask, maskBad, size, footprint);
     62    // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask
     63    int numCols = input->numCols, numRows = input->numRows; // Image dimensions
     64    psImage *stampMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // Mask to use for stamps
     65    for (int y = 0; y < numRows; y++) {
     66        for (int x = 0; x < numCols; x++) {
     67            stampMask->data.PS_TYPE_MASK_DATA[y][x] =
     68                (refRO->mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) ? MASK_BAD : 0;
     69        }
     70    }
    5171
    5272#if 0
    … …  
    6686#endif
    6787
    68     pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order,
    69                                                                  widths, orders); // Kernel basis functions
     88    pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order, widths, orders,
     89                                                                 inner, binning); // Kernel basis functions
    7090    psArray *stamps = NULL;             // Stamps for matching PSF
    7191    psVector *solution = NULL;          // Solution to match PSF
    … …  
    7393    int numRejected = -1;               // Number of rejected stamps in each iteration
    7494    for (int i = 0; i < iter && numRejected != 0; i++) {
    75         stamps = pmSubtractionFindStamps(stamps, refRO->image, subMask, threshold, spacing);
     95        stamps = pmSubtractionFindStamps(stamps, refRO->image, stampMask, MASK_BAD, MASK_STAMP,
     96                                         threshold, spacing, size + footprint);
    7697        if (!stamps) {
    7798            psError(PS_ERR_UNKNOWN, false, "Unable to find stamps on reference image.");
    … …  
    91112        }
    92113
    93         numRejected = pmSubtractionRejectStamps(stamps, refRO->image, inRO->image, subMask,
     114        numRejected = pmSubtractionRejectStamps(stamps, refRO->image, inRO->image, stampMask, MASK_STAMP,
    94115                                                solution, footprint, rej, kernels);
    95116        if (numRejected < 0) {
    … …  
    99120        psLogMsg("ppSub", PS_LOG_INFO, "%d stamps rejected on iteration %d.", numRejected, i);
    100121    }
    101 
    102     if (numRejected > 0) {
    103         solution = pmSubtractionSolveEquation(solution, stamps);
    104         if (!solution) {
    105             psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    106             goto ERROR;
    107         }
    108     }
     122    psFree(stampMask);
    109123
    110124    psImage *convImage = NULL, *convWeight = NULL, *convMask = NULL; // Convolved images
    111     if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask, refRO->image, refRO->weight, subMask,
    112                                maskBlank, solution, kernels)) {
     125    if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask,
     126                               refRO->image, refRO->weight, refRO->mask,
     127                               MASK_BAD, maskBlank, solution, kernels)) {
    113128        psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image.");
    114129        goto ERROR;
    115130    }
    116     psFree(subMask);
    117131
    118132    // Do the subtraction
    119133    if (reverse) {
    120         outRO->image = (psImage*)psBinaryOp(NULL, convImage, "-", inRO->image);
     134        (void)psBinaryOp(inRO->image, convImage, "-", inRO->image);
    121135    } else {
    122         outRO->image = (psImage*)psBinaryOp(NULL, inRO->image, "-", convImage);
     136        (void)psBinaryOp(inRO->image, inRO->image, "-", convImage);
    123137    }
    124     outRO->mask = (psImage*)psBinaryOp(NULL, convMask, "|", inRO->mask);
    125     if (convWeight) {
    126         outRO->weight = (psImage*)psBinaryOp(NULL, convWeight, "+", inRO->weight);
    127     }
     138    (void)psBinaryOp(inRO->mask, convMask, "|", inRO->mask);
     139    (void)psBinaryOp(inRO->weight, convWeight, "+", inRO->weight);
    128140
    129141    psFree(convImage);
    130142    psFree(convMask);
    131143    psFree(convWeight);
    132 
    133     outRO->data_exists = true;
    134     outCell->data_exists = true;
    135     outCell->parent->data_exists = true;
    136 
    137     if (!pmFPACopyConcepts(outCell->parent->parent, inRO->parent->parent->parent)) {
    138         psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
    139         return false;
    140     }
    141144
    142145#if 0
    … …  
    173176#endif
    174177
    175     psFree(kernels);
    176     psFree(stamps);
    177     psFree(solution);
    178178
    179     if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) {
    180         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
    181         pmFPACopy(photFile->fpa, inRO->parent->parent->parent);
     179        psFree(kernels);
     180        psFree(stamps);
     181        psFree(solution);
     182        return true;
    182183
    183         psphotReadout(config, view);
    184 
    185         pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
    186     }
    187 
    188     return true;
    189 
    190 ERROR:
    191     psFree(kernels);
    192     psFree(stamps);
    193     psFree(solution);
    194     return false;
     184    ERROR:
     185        psFree(kernels);
     186        psFree(stamps);
     187        psFree(solution);
     188        return false;
    195189}
Note: See TracChangeset for help on using the changeset viewer.