IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 29, 2009, 6:15:31 PM (17 years ago)
Author:
beaumont
Message:

merged with head

Location:
branches/cnb_branches/cnb_branch_20090301
Files:
21 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/cnb_branches/cnb_branch_20090301

  • branches/cnb_branches/cnb_branch_20090301/ppStack

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/Makefile.am

    r23352 r23594  
    2929        ppStackPrepare.c        \
    3030        ppStackConvolve.c       \
     31        ppStackCombinePrepare.c \
    3132        ppStackCombineInitial.c \
    3233        ppStackReject.c         \

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStack.h

    r23352 r23594  
    77#include <pslib.h>
    88#include <psmodules.h>
     9
     10#include "ppStackOptions.h"
    911
    1012// Mask values for inputs
    … …  
    106108/// Convolve image to match specified seeing
    107109bool ppStackMatch(pmReadout *readout,   // Readout to be convolved; replaced with output
    108                   psArray **regions,    // Array of regions used in each PSF matching, returned
    109                   psArray **kernels,    // Array of kernels used in each PSF matching, returned
    110                   float *chi2,          // Chi^2 from the stamps
    111                   float *weighting,     // Stack weighting (1/noise^2)
    112                   psArray *sources,     // Array of sources
    113                   const pmPSF *psf,     // Target PSF
    114                   psRandom *rng,        // Random number generator
     110                  ppStackOptions *options, // Options for stacking
     111                  int index,            // Index of image to match
    115112                  const pmConfig *config // Configuration
    116113    );
    … …  
    120117///
    121118/// Corrects the source PSF photometry to a common system.  Return the sum of the exposure times.
    122 float ppStackSourcesTransparency(const psArray *sourceLists, // Sources for each input
    123                                  psVector *inputMask, // Indicates bad input
    124                                  const pmFPAview *view, // View to readout
    125                                  const pmConfig *config // Configuration
     119bool ppStackSourcesTransparency(ppStackOptions *options, // Stacking options
     120                                const pmFPAview *view, // View to readout
     121                                const pmConfig *config // Configuration
    126122    );
    127123

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackArguments.c

    r23352 r23594  
    2020    fprintf(stderr, "\nPan-STARRS Image combination\n\n");
    2121    fprintf(stderr,
    22             "Usage: %s INPUTS.mdc OUTPUT_ROOT [-sources STAMPS.cmf | -stamps STAMPS.dat]\n"
     22            "Usage: %s -input INPUTS.mdc OUTPUT_ROOT [-sources STAMPS.cmf | -stamps STAMPS.dat]\n"
    2323            "where INPUTS.mdc contains various METADATAs, each with:\n"
    2424            "\tIMAGE(STR):     Image filename\n"
    … …  
    2626            "\tVARIANCE(STR):  Variance map filename\n"
    2727            "\tPSF(STR):       PSF filename\n"
    28             "\tSOURCES(STR):   Sources filename\n"
    29             "\tWEIGHTING(F32): Relative weighting to be applied\n",
     28            "\tSOURCES(STR):   Sources filename\n",
    3029            program);
    3130    fprintf(stderr, "\n");
    … …  
    137136        psArgumentRemove(argNum, &argc, argv);
    138137    }
    139 
    140138
    141139    psMetadata *arguments = config->arguments; // Command-line arguments
    … …  
    187185    psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "visualisation", false);
    188186
    189     if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) {
     187    if (argc == 1) {
    190188        usage(argv[0], arguments, config);
    191189    }
     190
     191    if ((argNum = psArgumentGet(argc, argv, "-input"))) {
     192        psArgumentRemove(argNum, &argc, argv);
     193        if (argNum >= argc) {
     194            usage(argv[0], arguments, config);
     195        }
     196
     197        unsigned int numBad = 0;                     // Number of bad lines
     198        psMetadata *inputs = psMetadataConfigRead(NULL, &numBad, argv[argNum], false); // Input file info
     199        if (!inputs || numBad > 0) {
     200            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to cleanly read MDC file with inputs.");
     201            return false;
     202        }
     203        psMetadataAddMetadata(arguments, PS_LIST_TAIL, "INPUTS", 0, "Metadata with input details", inputs);
     204        psFree(inputs);
     205
     206        psArgumentRemove(argNum, &argc, argv);
     207    }
     208
     209    if (!psArgumentParse(arguments, &argc, argv) || argc != 2) {
     210        usage(argv[0], arguments, config);
     211    }
     212
     213    psMetadataAddStr(arguments, PS_LIST_TAIL, "OUTPUT", 0, "Root name of the output image list", argv[1]);
    192214
    193215    const char *stampsName = psMetadataLookupStr(NULL, arguments, "-stamps"); // Name of stamps file
    194216    psMetadataAddStr(arguments, PS_LIST_TAIL, "STAMPS", 0, "Stamps file", stampsName);
    195 
    196     unsigned int numBad = 0;                     // Number of bad lines
    197     psMetadata *inputs = psMetadataConfigRead(NULL, &numBad, argv[1], false); // Information about inputs
    198     if (!inputs || numBad > 0) {
    199         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to cleanly read MDC file with inputs.");
    200         return false;
    201     }
    202     psMetadataAddMetadata(arguments, PS_LIST_TAIL, "INPUTS", 0, "Metadata with input details", inputs);
    203     psFree(inputs);
    204     psMetadataAddStr(arguments, PS_LIST_TAIL, "OUTPUT", 0, "Root name of the output image list", argv[2]);
    205217
    206218    valueArgStr(arguments, "-stats", "STATS", arguments);
    … …  
    295307    psTrace("ppStack", 1, "Done parsing arguments\n");
    296308
    297     // Dump configuration, now that's it's settled
    298     bool status;
    299     psString dump_file =  psMetadataLookupStr(&status, config->arguments, "DUMP_CONFIG");
    300     if (dump_file) {
    301         pmConfigCamerasCull(config, NULL);
    302         pmConfigRecipesCull(config, "PPSTACK,PPSUB,PPSTATS,PSPHOT,MASKS,JPEG");
    303 
    304         pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Input file
    305         pmConfigDump(config, input->fpa, dump_file);
    306     }
     309    pmConfigCamerasCull(config, NULL);
     310    pmConfigRecipesCull(config, "PPSTACK,PPSUB,PPSTATS,PSPHOT,MASKS,JPEG");
    307311
    308312    return true;

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackCamera.c

    r21366 r23594  
    1111#include "ppStack.h"
    1212
    13 
    14 #if 0
    15 // Define an output convolved image file
    16 static pmFPAfile *defineOutputConvolved(const char *name, // FPA file name
    17                                         pmFPA *fpa, // FPA to bind
    18                                         const pmConfig *config, // Configuration
    19                                         pmFPAfileType type // Expected type
    20     )
     13// Define a file
     14static pmFPAfile *defineFile(pmConfig *config, // Configuration
     15                             pmFPAfile *bind, // File to which to bind
     16                             const char *name, // Name of file rule
     17                             const char *filename, // Name of file
     18                             pmFPAfileType type // Type of file
     19                             )
    2120{
    22     pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, name);
    23     if (!file) {
    24         psError(PS_ERR_UNKNOWN, false, "Unable to define output convolved file %s", name);
     21
     22    psArray *dummy = psArrayAlloc(1);   // Dummy array of filenames for this FPA
     23    dummy->data[0] = psStringCopy(filename);
     24    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "FILENAMES", PS_META_REPLACE,
     25                       "Filenames for file rule definition", dummy);
     26    psFree(dummy);
     27
     28    bool found = false;             // Found the file?
     29    pmFPAfile *file = bind ? pmFPAfileBindFromArgs(&found, bind, config, name, "FILENAMES") :
     30        pmFPAfileDefineFromArgs(&found, config, name, "FILENAMES");
     31    if (!file || !found) {
     32        psError(PS_ERR_UNKNOWN, false, "Unable to define file %s from %s", name, filename);
    2533        return NULL;
    2634    }
    27     if (file->type != PM_FPA_FILE_IMAGE) {
    28         psError(PS_ERR_IO, true, "PPSTACK.OUTCONV is not of type %s", pmFPAfileStringFromType(type));
     35    if (file->type != type) {
     36        psError(PS_ERR_IO, true, "%s is not of type %s", name, pmFPAfileStringFromType(type));
    2937        return NULL;
    3038    }
    … …  
    3341}
    3442
    35 // Define an input convolved image file, using the output as a basis
    36 static pmFPAfile *defineInputConvolved(const char *inputName, // Input FPA file name
    37                                        pmFPAfile *outFile, // Corresponding output FPA file
    38                                        pmConfig *config, // Configuration
    39                                        pmFPAfileType type // Expected type
    40     )
    41 {
    42     pmFPAview *view = pmFPAviewAlloc(0); // View into sky cells
    43     view->chip = view->cell = view->readout = 0;
    44 
    45     psString imageName = pmFPANameFromRule(outFile->filerule, outFile->fpa, view);
    46     psArray *imageNames = psArrayAlloc(1);
    47     imageNames->data[0] = imageName;
    48     psMetadataAddArray(config->arguments, PS_LIST_TAIL, "INCONV.FILENAMES", PS_META_REPLACE,
    49                        "Filenames of input convolved image files", imageNames);
    50     psFree(imageNames);
    51     bool found = false;                 // Found the file?
    52     pmFPAfile *imageFile = pmFPAfileDefineFromArgs(&found, config, "PPSTACK.INCONV",
    53                                                    "INCONV.FILENAMES");
    54     psMetadataRemoveKey(config->arguments, "INCONV.FILENAMES");
    55     if (!imageFile || !found) {
    56         psError(PS_ERR_UNKNOWN, false, "Unable to define %s file", inputName);
    57         return NULL;
    58     }
    59     if (imageFile->type != type) {
    60         psError(PS_ERR_IO, true, "PPSTACK.INCONV is not of type %s",
    61                 pmFPAfileStringFromType(type));
    62         return NULL;
    63     }
    64 
    65     return imageFile;
    66 }
    67 #endif
    68 
    6943
    7044
    7145bool ppStackCamera(pmConfig *config)
    7246{
    73     bool haveVariances = false;           // Do we have variance maps?
     47    int num = 0;                        // Number of inputs
     48    bool haveVariances = false;         // Do we have variance maps?
    7449    bool havePSFs = false;              // Do we have PSFs?
    7550
    76     psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info
    77     psMetadataIterator *iter = psMetadataIteratorAlloc(inputs, PS_LIST_HEAD, NULL); // Iterator
    78     psMetadataItem *item;               // Item from iteration
    79     int i = 0;                          // Counter
    80     while ((item = psMetadataGetAndIncrement(iter))) {
    81         if (item->type != PS_DATA_METADATA) {
    82             psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    83                     "Component %s of the input metadata is not of type METADATA", item->name);
    84             psFree(iter);
    85             return false;
    86         }
    87 
    88         psMetadata *input = item->data.md; // The input metadata of interest
    89 
    90         psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image
    91         if (!image || strlen(image) == 0) {
    92             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name);
    93             psFree(iter);
    94             return false;
    95         }
    96 
    97         bool mdok;
    98         psString mask = psMetadataLookupStr(&mdok, input, "MASK"); // Name of mask
    99         psString variance = psMetadataLookupStr(&mdok, input, "VARIANCE"); // Name of variance map
    100         psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF
    101         psString sources = psMetadataLookupStr(&mdok, input, "SOURCES"); // Name of sources
    102 
    103         // Add the image file
    104         psArray *imageFiles = psArrayAlloc(1); // Array of filenames for this FPA
    105         imageFiles->data[0] = psMemIncrRefCounter(image);
    106         psMetadataAddArray(config->arguments, PS_LIST_TAIL, "IMAGE.FILENAMES", PS_META_REPLACE,
    107                            "Filenames of image files", imageFiles);
    108         psFree(imageFiles);
    109 
    110         bool found = false;             // Found the file?
    111         pmFPAfile *imageFile = pmFPAfileDefineFromArgs(&found, config, "PPSTACK.INPUT", "IMAGE.FILENAMES");
    112         if (!imageFile || !found) {
    113             psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image);
    114             return false;
    115         }
    116         if (imageFile->type != PM_FPA_FILE_IMAGE) {
    117             psError(PS_ERR_IO, true, "PPSTACK.INPUT is not of type IMAGE");
    118             return false;
    119         }
    120 
    121         // Optionally add the mask file
    122         if (mask && strlen(mask) > 0) {
    123             psArray *maskFiles = psArrayAlloc(1); // Array of filenames for this FPA
    124             maskFiles->data[0] = psMemIncrRefCounter(mask);
    125             psMetadataAddArray(config->arguments, PS_LIST_TAIL, "MASK.FILENAMES", PS_META_REPLACE,
    126                                "Filenames of mask files", maskFiles);
    127             psFree(maskFiles);
    128 
    129             bool status;
    130             pmFPAfile *maskFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPSTACK.INPUT.MASK",
    131                                                         "MASK.FILENAMES");
     51    bool status = false;                // Status of file definition
     52
     53    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppSim
     54    if (!recipe) {
     55        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE);
     56        return false;
     57    }
     58    bool convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images before stack?
     59
     60    psArray *runImages = pmFPAfileDefineMultipleFromRun(&status, NULL, config,
     61                                                        "PPSTACK.INPUT"); // Input images from previous run
     62    if (runImages) {
     63        // Defining files from the RUN metadata
     64        num = runImages->n;
     65
     66        psArray *runMasks = pmFPAfileDefineMultipleFromRun(&status, runImages, config,
     67                                                           "PPSTACK.INPUT.MASK"); // Input masks
     68        if (!status) {
     69            psError(PS_ERR_UNKNOWN, false, "Unable to define input masks from RUN metadata.");
     70            psFree(runImages);
     71            return false;
     72        }
     73        psFree(runMasks);
     74
     75        psArray *runVars = pmFPAfileDefineMultipleFromRun(&status, runImages, config,
     76                                                          "PPSTACK.INPUT.VARIANCE"); // Input variances
     77        if (!status) {
     78            psError(PS_ERR_UNKNOWN, false, "Unable to define input variances from RUN metadata.");
     79            psFree(runImages);
     80            return false;
     81        }
     82        if (runVars) {
     83            haveVariances = true;
     84        }
     85        psFree(runVars);
     86
     87        psArray *runPSF = pmFPAfileDefineMultipleFromRun(&status, runImages, config,
     88                                                         "PPSTACK.INPUT.PSF"); // Input PSFs
     89        if (!status) {
     90            psError(PS_ERR_UNKNOWN, false, "Unable to define input PSFs from RUN metadata.");
     91            psFree(runImages);
     92            return false;
     93        }
     94        if (runPSF) {
     95            havePSFs = true;
     96        }
     97        psFree(runPSF);
     98
     99        psArray *runSrc = pmFPAfileDefineMultipleFromRun(&status, runImages, config,
     100                                                         "PPSTACK.INPUT.SOURCES"); // Input sources
     101        if (!status) {
     102            psError(PS_ERR_UNKNOWN, false, "Unable to define input sources from RUN metadata.");
     103            psFree(runImages);
     104            return false;
     105        }
     106        if (!runSrc) {
     107            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to define input sources from RUN metadata.");
     108            psFree(runImages);
     109            return false;
     110        }
     111        psFree(runSrc);
     112
     113        if (convolve) {
     114            psArray *runKernel = pmFPAfileDefineMultipleFromRun(&status, runImages, config,
     115                                                                "PPSTACK.CONV.KERNEL"); // Convolution kernels
    132116            if (!status) {
     117                psError(PS_ERR_UNKNOWN, false, "Unable to define convolution kernels from RUN metadata.");
     118                psFree(runImages);
     119                return false;
     120            }
     121            if (!runKernel) {
     122                psError(PS_ERR_UNEXPECTED_NULL, true,
     123                        "Unable to define convolution kernels from RUN metadata.");
     124                psFree(runImages);
     125                return false;
     126            }
     127            psFree(runKernel);
     128        }
     129
     130        psFree(runImages);
     131    } else {
     132        // Defining files from the input metadata
     133        psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info
     134        if (!inputs) {
     135            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find inputs.");
     136            return false;
     137        }
     138        psMetadataIterator *iter = psMetadataIteratorAlloc(inputs, PS_LIST_HEAD, NULL); // Iterator
     139        psMetadataItem *item;               // Item from iteration
     140        int i = 0;                          // Counter
     141        while ((item = psMetadataGetAndIncrement(iter))) {
     142            if (item->type != PS_DATA_METADATA) {
     143                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
     144                        "Component %s of the input metadata is not of type METADATA", item->name);
     145                psFree(iter);
     146                return false;
     147            }
     148
     149            psMetadata *input = item->data.md; // The input metadata of interest
     150
     151            psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image
     152            if (!image || strlen(image) == 0) {
     153                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name);
     154                psFree(iter);
     155                return false;
     156            }
     157
     158            bool mdok;
     159            psString mask = psMetadataLookupStr(&mdok, input, "MASK"); // Name of mask
     160            psString variance = psMetadataLookupStr(&mdok, input, "VARIANCE"); // Name of variance map
     161            psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF
     162            psString sources = psMetadataLookupStr(&mdok, input, "SOURCES"); // Name of sources
     163
     164            pmFPAfile *imageFile = defineFile(config, NULL, "PPSTACK.INPUT",
     165                                              image, PM_FPA_FILE_IMAGE); // File for image
     166            if (!imageFile) {
     167                psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image);
     168                return false;
     169            }
     170
     171            if (mask && strlen(mask) > 0 &&
     172                !defineFile(config, imageFile, "PPSTACK.INPUT.MASK", mask, PM_FPA_FILE_MASK)) {
    133173                psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
    134174                return false;
    135175            }
    136             if (maskFile->type != PM_FPA_FILE_MASK) {
    137                 psError(PS_ERR_IO, true, "PPSTACK.INPUT.MASK is not of type MASK");
    138                 return false;
    139             }
    140         }
    141 
    142         // Optionally add the variance file
    143         if (variance && strlen(variance) > 0) {
    144             haveVariances = true;
    145             psArray *varianceFiles = psArrayAlloc(1); // Array of filenames for this FPA
    146             varianceFiles->data[0] = psMemIncrRefCounter(variance);
    147             psMetadataAddArray(config->arguments, PS_LIST_TAIL, "VARIANCE.FILENAMES", PS_META_REPLACE,
    148                                "Filenames of variance files", varianceFiles);
    149             psFree(varianceFiles);
    150 
    151             bool status;
    152             pmFPAfile *varianceFile = pmFPAfileBindFromArgs(&status, imageFile, config,
    153                                                             "PPSTACK.INPUT.VARIANCE", "VARIANCE.FILENAMES");
    154             if (!status) {
    155                 psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance);
    156                 return false;
    157             }
    158             if (varianceFile->type != PM_FPA_FILE_VARIANCE) {
    159                 psError(PS_ERR_IO, true, "PPSTACK.INPUT.VARIANCE is not of type VARIANCE");
    160                 return false;
    161             }
    162         }
    163 
    164         // Add the psf file
    165         if (!psf || strlen(psf) == 0) {
    166             if (havePSFs) {
    167                 psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF %d", i);
    168                 return false;
    169             }
    170         } else {
    171             if (!havePSFs && i != 0) {
    172                 psWarning("PSF not provided for all inputs --- ignoring.");
    173             } else {
    174                 psArray *psfFiles = psArrayAlloc(1); // Array of filenames for this FPA
    175                 psfFiles->data[0] = psMemIncrRefCounter(psf);
    176                 psMetadataAddArray(config->arguments, PS_LIST_TAIL, "PSF.FILENAMES", PS_META_REPLACE,
    177                                    "Filenames of PSF files", psfFiles);
    178                 psFree(psfFiles);
    179 
    180                 bool status;
    181                 pmFPAfile *psfFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPSTACK.INPUT.PSF",
    182                                                            "PSF.FILENAMES");
    183                 if (!status) {
    184                     psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf);
     176
     177            if (variance && strlen(variance) > 0) {
     178                haveVariances = true;
     179                if (!defineFile(config, imageFile, "PPSTACK.INPUT.VARIANCE", variance,
     180                                PM_FPA_FILE_VARIANCE)) {
     181                    psError(PS_ERR_UNKNOWN, false,
     182                            "Unable to define file from variance %d (%s)", i, variance);
    185183                    return false;
    186184                }
    187                 if (psfFile->type != PM_FPA_FILE_PSF) {
    188                     psError(PS_ERR_IO, true, "PPSTACK.INPUT.PSF is not of type PSF");
     185            }
     186
     187            if (psf && strlen(psf) > 0) {
     188                if (i != 0 && !havePSFs) {
     189                    psWarning("PSF not provided for all inputs --- ignoring.");
     190                } else {
     191                    havePSFs = true;
     192                    if (!defineFile(config, imageFile, "PPSTACK.INPUT.PSF", psf, PM_FPA_FILE_PSF)) {
     193                        psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf);
     194                        return false;
     195                    }
     196                }
     197            } else if (havePSFs) {
     198                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF %d", i);
     199                return false;
     200            }
     201
     202            if (!sources || strlen(sources) == 0) {
     203                psError(PS_ERR_UNEXPECTED_NULL, true, "SOURCES not provided for file %d", i);
     204                return false;
     205            }
     206            if (!defineFile(config, imageFile, "PPSTACK.INPUT.SOURCES", sources, PM_FPA_FILE_CMF)) {
     207                psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)",
     208                        i, sources);
     209                return false;
     210            }
     211
     212            if (convolve) {
     213                pmFPAfile *kernel = pmFPAfileDefineOutput(config, imageFile->fpa, "PPSTACK.CONV.KERNEL");
     214                if (!kernel) {
     215                    psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.CONV.KERNEL"));
    189216                    return false;
    190217                }
    191                 havePSFs = true;
    192             }
    193         }
    194 
    195         // Add the sources file
    196         {
    197             psArray *sourcesFiles = psArrayAlloc(1); // Array of filenames for this FPA
    198             sourcesFiles->data[0] = psMemIncrRefCounter(sources);
    199             psMetadataAddArray(config->arguments, PS_LIST_TAIL, "SOURCES.FILENAMES", PS_META_REPLACE,
    200                                "Filenames of SOURCES files", sourcesFiles);
    201             psFree(sourcesFiles);
    202 
    203             bool status;
    204             pmFPAfile *sourcesFile = pmFPAfileBindFromArgs(&status, imageFile, config,
    205                                                            "PPSTACK.INPUT.SOURCES", "SOURCES.FILENAMES");
    206             if (!status) {
    207                 psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)",
    208                         i, sources);
    209                 return false;
    210             }
    211             if (sourcesFile->type != PM_FPA_FILE_CMF) {
    212                 psError(PS_ERR_IO, true, "PPSTACK.INPUT.SOURCES is not of type CMF");
    213                 return false;
    214             }
    215         }
    216 
    217 #if 0
    218         // Output convolved files
    219         pmFPAfile *outconvImage    = defineOutputConvolved("PPSTACK.OUTCONV", imageFile->fpa, config,
    220                                                            PM_FPA_FILE_IMAGE);
    221         pmFPAfile *outconvMask     = defineOutputConvolved("PPSTACK.OUTCONV.MASK", imageFile->fpa, config,
    222                                                            PM_FPA_FILE_MASK);
    223         pmFPAfile *outconvVariance = defineOutputConvolved("PPSTACK.OUTCONV.VARIANCE", imageFile->fpa, config,
    224                                                            PM_FPA_FILE_VARIANCE);
    225         if (!outconvImage || !outconvMask || !outconvVariance) {
    226             return false;
    227         }
    228 
    229         // Input convolved files
    230         pmFPAfile *inconvImage    = defineInputConvolved("PPSTACK.INCONV", outconvImage, config,
    231                                                          PM_FPA_FILE_IMAGE);
    232         pmFPAfile *inconvMask     = defineInputConvolved("PPSTACK.INCONV.MASK", outconvMask, config,
    233                                                          PM_FPA_FILE_MASK);
    234         pmFPAfile *inconvVariance = defineInputConvolved("PPSTACK.INCONV.VARIANCE", outconvVariance, config,
    235                                                          PM_FPA_FILE_VARIANCE);
    236         if (!inconvImage || !inconvMask || !inconvVariance) {
    237             return false;
    238         }
    239 #endif
    240 
    241         i++;
    242     }
    243     psFree(iter);
    244     psMetadataRemoveKey(config->arguments, "IMAGE.FILENAMES");
    245     if (psMetadataLookup(config->arguments, "MASK.FILENAMES")) {
    246         psMetadataRemoveKey(config->arguments, "MASK.FILENAMES");
    247     }
    248     if (psMetadataLookup(config->arguments, "VARIANCE.FILENAMES")) {
    249         psMetadataRemoveKey(config->arguments, "VARIANCE.FILENAMES");
    250     }
    251     if (psMetadataLookup(config->arguments, "PSF.FILENAMES")) {
    252         psMetadataRemoveKey(config->arguments, "PSF.FILENAMES");
    253     }
    254     if (psMetadataLookup(config->arguments, "SOURCES.FILENAMES")) {
    255         psMetadataRemoveKey(config->arguments, "SOURCES.FILENAMES");
    256     }
    257 
    258     psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.NUM", 0, "Number of input files", i);
     218                kernel->save = true;
     219            }
     220
     221            i++;
     222        }
     223        psFree(iter);
     224        psMetadataRemoveKey(config->arguments, "FILENAMES");
     225        num = i;
     226    }
     227    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.NUM", 0, "Number of input files", num);
    259228    psMetadataAddBool(config->arguments, PS_LIST_TAIL, "HAVE.PSF", 0, "Have PSFs available?", havePSFs);
    260229
    … …  
    343312    }
    344313    jpeg2->save = true;
    345 
    346 
    347     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppSim
    348     if (!recipe) {
    349         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE);
    350         return false;
    351     }
    352314
    353315    // For photometry, we operate on the chip-mosaicked image

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackCombineInitial.c

    r23352 r23594  
    1616    psAssert(config, "Require configuration");
    1717
    18     psTimerStart("PPSTACK_FINAL");
     18    if (!options->convolve) {
     19        // No need to do initial combination when we haven't convolved
     20        return true;
     21    }
    1922
    20     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    21     psAssert(recipe, "We've thrown an error on this before.");
     23    psTimerStart("PPSTACK_INITIAL");
    2224
    2325    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    … …  
    2527                                          "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans
    2628
    27     pmFPAview *view = NULL;             // View to readout
    28 
    29     int row0, col0;                 // Offset for readout
    30     int numCols, numRows;           // Size of readout
    31     ppStackThread *thread = stack->threads->data[0]; // Representative thread
    32     if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) {
    33         psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
    34         return false;
    35     }
    36 
    37     pmFPAfileActivate(config->files, false, NULL);
    38     ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true);
    39     view = ppStackFilesIterateDown(config);
    40     if (!view) {
    41         return false;
    42     }
    43 
    44     pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
    45     options->outRO = pmReadoutAlloc(outCell); // Output readout
    46     psFree(view);
    47 
    48     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    49     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    50     if (!pmReadoutStackDefineOutput(options->outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
    51         psError(PS_ERR_UNKNOWN, false, "Unable to prepare output.");
    52         return false;
    53     }
    54 
    55     psFree(options->cells); options->cells = NULL;
    5629
    5730    bool status;                    // Status of read
    … …  
    6942        }
    7043
    71         // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)
    7244        psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start
    7345        psArrayAdd(job->args, 1, thread);
    … …  
    128100
    129101    if (options->stats) {
    130         psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_FINAL", 0,
    131                          "Time to make final stack", psTimerMark("PPSTACK_FINAL"));
     102        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_INITIAL", 0,
     103                         "Time to make final stack", psTimerMark("PPSTACK_INITIAL"));
    132104    }
    133105

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackConvolve.c

    r23352 r23594  
    3030    int num = options->num;             // Number of inputs
    3131    options->cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again
    32     options->subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking
    33     options->subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking
     32    options->kernels = psArrayAlloc(num); // PSF-matching kernels --- required in the stacking
     33    options->regions = psArrayAlloc(num); // PSF-matching regions --- required in the stacking
    3434    int numGood = 0;                    // Number of good frames
    3535    options->numCols = 0;
    … …  
    3838    psVectorInit(options->matchChi2, NAN);
    3939    options->weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
    40     psVectorInit(options->weightings, NAN);
     40    psVectorInit(options->weightings, 0.0);
    4141    options->covariances = psArrayAlloc(num); // Covariance matrices
    4242
    … …  
    7878
    7979        // Background subtraction, scaling and normalisation is performed automatically by the image matching
    80         psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
    8180        psTimerStart("PPSTACK_MATCH");
    82 
    83         if (!ppStackMatch(readout, &regions, &kernels, &options->matchChi2->data.F32[i],
    84                           &options->weightings->data.F32[i], options->sourceLists->data[i],
    85                           options->psf, rng, config)) {
     81        if (!ppStackMatch(readout, options, i, config)) {
    8682            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
    87             options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_MATCH;
     83            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH;
    8884            psErrorClear();
    8985            continue;
    … …  
    9288
    9389        if (options->stats) {
    94             pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
    95                                                                 PM_SUBTRACTION_ANALYSIS_KERNEL);// Conv kernel
    9690            psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_MATCH", PS_META_DUPLICATE_OK,
    9791                             "Time to match PSF", psTimerMark("PPSTACK_MATCH"));
    98             psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK,
    99                              "Mean deviation for stamps", kernels->mean);
    100             psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK,
    101                              "RMS deviation for stamps", kernels->rms);
    102             psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK,
    103                              "Number of stamps", kernels->numStamps);
    104             float deconv = psMetadataLookupF32(NULL, readout->analysis,
    105                                                PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Deconvolution fraction
    106             psMetadataAddF32(options->stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK,
    107                              "Deconvolution fraction for kernel", deconv);
    10892            psMetadataAddF32(options->stats, PS_LIST_TAIL, "PPSTACK.WEIGHTING", PS_META_DUPLICATE_OK,
    10993                             "Weighting for image", options->weightings->data.F32[i]);
     94
     95            if (options->convolve) {
     96                // Pull parameters out of convolution kernel
     97                pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
     98                                                                    PM_SUBTRACTION_ANALYSIS_KERNEL);
     99                psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK,
     100                                 "Mean deviation for stamps", kernels->mean);
     101                psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK,
     102                                 "RMS deviation for stamps", kernels->rms);
     103                psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK,
     104                                 "Number of stamps", kernels->numStamps);
     105                float deconv = psMetadataLookupF32(NULL, readout->analysis,
     106                                                   PM_SUBTRACTION_ANALYSIS_DECONV_MAX);
     107                psMetadataAddF32(options->stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK,
     108                                 "Deconvolution fraction for kernel", deconv);
     109            }
    110110        }
    111111        psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH"));
    112112
    113         options->subRegions->data[i] = regions;
    114         options->subKernels->data[i] = kernels;
    115 
    116113        // Write the temporary convolved files
    117         pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
    118         assert(hdu);
    119         ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config);
    120         psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
    121         pmConfigMaskWriteHeader(config, maskHeader);
    122         ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config);
    123         psFree(maskHeader);
    124         psImageCovarianceTransfer(readout->variance, readout->covariance);
    125         ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config);
     114        if (options->convolve) {
     115            pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
     116            assert(hdu);
     117            ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config);
     118            psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
     119            pmConfigMaskWriteHeader(config, maskHeader);
     120            ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config);
     121            psFree(maskHeader);
     122            psImageCovarianceTransfer(readout->variance, readout->covariance);
     123            ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config);
    126124#ifdef TESTING
    127         {
    128             psString name = NULL;
    129             psStringAppend(&name, "covariance_%d.fits", i);
    130             ppStackWriteImage(name, hdu->header, readout->covariance->image, config);
    131             pmStackVisualPlotTestImage(readout->covariance->image, name);
    132             psFree(name);
    133         }
     125            {
     126                psString name = NULL;
     127                psStringAppend(&name, "covariance_%d.fits", i);
     128                ppStackWriteImage(name, hdu->header, readout->covariance->image, config);
     129                pmStackVisualPlotTestImage(readout->covariance->image, name);
     130                psFree(name);
     131            }
    134132#endif
     133        }
    135134
    136135        pmCell *inCell = readout->parent; // Input cell
    … …  
    152151    psFree(rng);
    153152
     153    psFree(options->norm); options->norm = NULL;
    154154    psFree(options->sourceLists); options->sourceLists = NULL;
    155155    psFree(options->psf); options->psf = NULL;
    … …  
    183183
    184184    // Reject images out-of-hand on the basis of their match chi^2
    185     {
     185    if (options->convolve) {
    186186        psVector *values = psVectorAllocEmpty(num, PS_TYPE_F32); // Values to sort
    187187        for (int i = 0; i < num; i++) {

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackFiles.c

    r23352 r23594  
    1717
    1818/// Files required for the convolution
    19 static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL };
     19static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE",
     20                                 "PPSTACK.CONV.KERNEL", NULL };
    2021
    2122/// Output files for the combination

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackFinish.c

    r23352 r23594  
    6060    }
    6161
     62
     63    // Dump configuration
     64    bool mdok;                          // Status of MD lookup
     65    psString dump = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG"); // File for config
     66    if (dump) {
     67        pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Input file
     68        pmConfigDump(config, input->fpa, dump);
     69    }
     70
    6271    return true;
    6372}

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackLoop.c

    r23352 r23594  
    6262        return false;
    6363    }
     64    psFree(options->cells); options->cells = NULL;
     65
     66    // Prepare for combination
     67    if (!ppStackCombinePrepare(stack, options, config)) {
     68        psError(PS_ERR_UNKNOWN, false, "Unable to prepare for combination.");
     69        psFree(stack);
     70        psFree(options);
     71        return false;
     72    }
    6473
    6574    // Initial combination
    … …  
    103112    // Clean up
    104113    psTrace("ppStack", 2, "Cleaning up after combination....\n");
    105     if (!ppStackCombineFinal(stack, options, config)) {
     114    if (!ppStackCleanup(stack, options, config)) {
    106115        psError(PS_ERR_UNKNOWN, false, "Unable to clean up.");
    107116        psFree(stack);

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackLoop.h

    r23352 r23594  
    3434    );
    3535
     36// Prepare for combination
     37bool ppStackCombinePrepare(
     38    ppStackThreadData *stack,           // Stack
     39    ppStackOptions *options,            // Options
     40    pmConfig *config                    // Configuration
     41    );
     42
    3643// Initial combination
    3744bool ppStackCombineInitial(
    … …  
    5461    );
    5562
     63// Cleanup following combination
     64bool ppStackCleanup(
     65    ppStackThreadData *stack,           // Stack
     66    ppStackOptions *options,            // Options
     67    pmConfig *config                    // Configuration
     68    );
     69
    5670// Photometry
    5771bool ppStackPhotometry(

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackMatch.c

    r23352 r23594  
    163163
    164164
    165 bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
    166                   psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
     165bool ppStackMatch(pmReadout *readout, ppStackOptions *options, int index, const pmConfig *config)
    167166{
    168167    assert(readout);
    169     assert(regions && !*regions);
    170     assert(kernels && !*kernels);
     168    assert(options);
    171169    assert(config);
    172     *weighting = 0.0;
    173170
    174171    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    … …  
    197194    }
    198195
    199     pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
    200 
    201     static int numInput = -1;            // Index of input file
    202     numInput++;
     196    // Match the PSF
     197    if (options->convolve) {
     198        pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily
    203199#ifdef TESTING
    204     // Read previously produced kernel
    205     if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
    206         const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
    207         assert(outName);
    208         // Read convolution kernel
    209         psString filename = NULL;   // Output filename
    210         psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
    211         psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
    212         psFree(filename);
    213         psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
    214         psFree(resolved);
    215         if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
    216             psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
     200        // This is a hack to use the temporary convolved images and kernel generated previously.
     201        // This makes the 'matching' operation much faster, allowing debugging of the stack process easier.
     202        // It implicitly assumes the output root name is the same between invocations.
     203
     204        // Read previously produced kernel
     205        if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
     206            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
     207            psAssert(file, "Require file");
     208
     209            pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
     210            view->chip = view->cell = view->readout = 0;
     211            psString filename = pmFPAfileNameFromRule(filerule->rule, file, view); // Filename of interest
     212
     213            // Read convolution kernel
     214            psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
     215            psFree(filename);
     216            psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
     217            psFree(resolved);
     218            if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) {
     219                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
     220                psFitsClose(fits);
     221                return false;
     222            }
    217223            psFitsClose(fits);
    218             return false;
    219         }
    220         psFitsClose(fits);
    221 
    222         // Add in variance factor
    223         pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
    224                                                             PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
    225         float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
    226         psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
    227         if (!isfinite(vf)) {
    228             vf = 1.0;
    229         }
    230         if (isfinite(vfItem->data.F32)) {
    231             vfItem->data.F32 *= vf;
    232         } else {
    233             vfItem->data.F32 = vf;
    234         }
    235 
    236         // Read image, mask, variance
    237         const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for image
    238         const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for mask
    239         const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for variance map
    240         psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
    241         psStringAppend(&imageName, "%s.%d.%s", outName, numInput, tempImage);
    242         psStringAppend(&maskName, "%s.%d.%s", outName, numInput, tempMask);
    243         psStringAppend(&varianceName, "%s.%d.%s", outName, numInput, tempVariance);
    244 
    245         if (!readImage(&readout->image, imageName, config) || !readImage(&readout->mask, maskName, config) ||
    246             !readImage(&readout->variance, varianceName, config)) {
    247             psError(PS_ERR_IO, false, "Unable to read previously produced image.");
     224
     225            // Add in variance factor
     226            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis,
     227                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
     228            float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
     229            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
     230            if (!isfinite(vf)) {
     231                vf = 1.0;
     232            }
     233            if (isfinite(vfItem->data.F32)) {
     234                vfItem->data.F32 *= vf;
     235            } else {
     236                vfItem->data.F32 = vf;
     237            }
     238
     239            if (!readImage(&readout->image, options->imageNames->data[index], config) ||
     240                !readImage(&readout->mask, options->maskNames->data[index], config) ||
     241                !readImage(&readout->variance, options->varianceNames->data[index], config)) {
     242                psError(PS_ERR_IO, false, "Unable to read previously produced image.");
     243                psFree(imageName);
     244                psFree(maskName);
     245                psFree(varianceName);
     246                return false;
     247            }
    248248            psFree(imageName);
    249249            psFree(maskName);
    250250            psFree(varianceName);
    251             return false;
    252         }
    253         psFree(imageName);
    254         psFree(maskName);
    255         psFree(varianceName);
    256 
    257         psRegion *region = psMetadataLookupPtr(NULL, output->analysis,
    258                                                PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
    259 
    260         pmSubtractionAnalysis(readout->analysis, kernels, region,
    261                               readout->image->numCols, readout->image->numRows);
    262 
    263         psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
    264         psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // New covariance matrix
    265         psFree(readout->covariance);
    266         readout->covariance = covar;
    267         psFree(kernel);
    268 
    269     } else {
    270 #endif
    271 
    272         // Normal operations here
    273         if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) {
    274             assert(psf);
    275             assert(sources);
     251
     252            psRegion *region = psMetadataLookupPtr(NULL, conv->analysis,
     253                                                   PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
     254
     255            pmSubtractionAnalysis(readout->analysis, kernels, region,
     256                                  readout->image->numCols, readout->image->numRows);
     257
     258            psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
     259            psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
     260            psFree(readout->covariance);
     261            readout->covariance = covar;
     262            psFree(kernel);
     263        } else {
     264#endif
     265
     266            // Normal operations here
     267            psAssert(options->psf, "Require target PSF");
     268            psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list");
    276269
    277270            int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order
    … …  
    284277            float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
    285278            float sysError = psMetadataLookupF32(NULL, ppsub, "SYS"); // Relative systematic error in kernel
    286             pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(
    287                 psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE")); // Kernel type
     279            const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type
     280            pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
    288281            psVector *widths = psMetadataLookupPtr(NULL, ppsub, "ISIS.WIDTHS"); // ISIS Gaussian widths
    289282            psVector *orders = psMetadataLookupPtr(NULL, ppsub, "ISIS.ORDERS"); // ISIS Polynomial orders
    … …  
    308301            }
    309302
    310 #if 0
    311             // Testing the normalisation of the fake image
    312             {
    313                 pmReadout *test = pmReadoutAlloc(NULL); // Test readout
    314                 psArray *sources = psArrayAlloc(1);     // Array of sources
    315                 pmSource *source = pmSourceAlloc();     // Source
    316                 sources->data[0] = source;
    317                 source->peak = pmPeakAlloc(500, 500, 10000, PM_PEAK_LONE);
    318                 source->psfMag = -13.0;
    319                 pmReadoutFakeFromSources(test, 1000, 1000, sources, NULL, NULL, psf, 0.1, 0, false, true);
    320                 float sum = 0.0;
    321                 for (int y = 0; y < test->image->numRows; y++) {
    322                     for (int x = 0; x < test->image->numCols; x++) {
    323                         sum += test->image->data.F32[y][x];
    324                     }
    325                 }
    326                 fprintf(stderr, "Photometry: %f --> %f = -13.0 ???\n", sum, -2.5*log10(sum));
    327 
    328                 psFits *fits = psFitsOpen("testphot.fits", "w");
    329                 psFitsWriteImage(fits, NULL, test->image, 0, NULL);
    330                 psFitsClose(fits);
    331                 exit(0);
    332             }
    333 #endif
    334 
    335303            pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
    336304
    337305            // For the sake of stamps, remove nearby sources
    338             psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources
     306            psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index],
     307                                                       footprint); // Filtered list of sources
    339308
    340309            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
    341                                           stampSources, NULL, NULL, psf, NAN, footprint + size,
     310                                          stampSources, NULL, NULL, options->psf, NAN, footprint + size,
    342311                                          false, true)) {
    343312                psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
    344313                psFree(fake);
    345314                psFree(optWidths);
    346                 psFree(output);
     315                psFree(conv);
    347316                return false;
    348317            }
    … …  
    357326                pmHDU *hdu = pmHDUFromCell(readout->parent);
    358327                psString name = NULL;
    359                 psStringAppend(&name, "fake_%03d.fits", numInput);
     328                psStringAppend(&name, "fake_%03d.fits", index);
    360329                pmStackVisualPlotTestImage(fake->image, name);
    361330                psFits *fits = psFitsOpen(name, "w");
    … …  
    367336                pmHDU *hdu = pmHDUFromCell(readout->parent);
    368337                psString name = NULL;
    369                 psStringAppend(&name, "real_%03d.fits", numInput);
     338                psStringAppend(&name, "real_%03d.fits", index);
    370339                pmStackVisualPlotTestImage(readout->image, name);
    371340                psFits *fits = psFitsOpen(name, "w");
    … …  
    381350
    382351            // Do the image matching
    383             if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, stride, regionSize, spacing,
    384                                     threshold, stampSources, stampsName, type, size, order, widths, orders,
    385                                     inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder,
    386                                     optThresh, iter, rej, sysError, maskVal, maskBad, maskPoor, poorFrac,
    387                                     badFrac, PM_SUBTRACTION_MODE_1)) {
    388                 psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
    389                 psFree(fake);
    390                 psFree(optWidths);
    391                 psFree(stampSources);
    392                 psFree(output);
    393                 return false;
     352            pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readout->analysis,
     353                                                               PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
     354            if (kernel) {
     355                if (!pmSubtractionMatchPrecalc(conv, NULL, readout, fake, readout->analysis,
     356                                               stride, sysError, maskVal, maskBad, maskPoor,
     357                                               poorFrac, badFrac)) {
     358                    psError(PS_ERR_UNKNOWN, false, "Unable to convolve images.");
     359                    psFree(fake);
     360                    psFree(optWidths);
     361                    psFree(stampSources);
     362                    psFree(conv);
     363                    return false;
     364                }
     365            } else {
     366                if (!pmSubtractionMatch(conv, NULL, readout, fake, footprint, stride, regionSize, spacing,
     367                                        threshold, stampSources, stampsName, type, size, order, widths,
     368                                        orders, inner, ringsOrder, binning, penalty,
     369                                        optimum, optWidths, optOrder, optThresh, iter, rej, sysError,
     370                                        maskVal, maskBad, maskPoor, poorFrac, badFrac,
     371                                        PM_SUBTRACTION_MODE_1)) {
     372                    psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
     373                    psFree(fake);
     374                    psFree(optWidths);
     375                    psFree(stampSources);
     376                    psFree(conv);
     377                    return false;
     378                }
    394379            }
    395380
    … …  
    398383                pmHDU *hdu = pmHDUFromCell(readout->parent);
    399384                psString name = NULL;
    400                 psStringAppend(&name, "conv_%03d.fits", numInput);
    401                 pmStackVisualPlotTestImage(output->image, name);
     385                psStringAppend(&name, "conv_%03d.fits", index);
     386                pmStackVisualPlotTestImage(conv->image, name);
    402387                psFits *fits = psFitsOpen(name, "w");
    403388                psFree(name);
    404                 psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
     389                psFitsWriteImage(fits, hdu->header, conv->image, 0, NULL);
    405390                psFitsClose(fits);
    406391            }
    … …  
    408393                pmHDU *hdu = pmHDUFromCell(readout->parent);
    409394                psString name = NULL;
    410                 psStringAppend(&name, "diff_%03d.fits", numInput);
     395                psStringAppend(&name, "diff_%03d.fits", index);
    411396                pmStackVisualPlotTestImage(fake->image, name);
    412397                psFits *fits = psFitsOpen(name, "w");
    413398                psFree(name);
    414                 psBinaryOp(fake->image, output->image, "-", fake->image);
     399                psBinaryOp(fake->image, conv->image, "-", fake->image);
    415400                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
    416401                psFitsClose(fits);
    … …  
    428413            // Set the variance factor
    429414            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
    430             float vf = psMetadataLookupF32(NULL, output->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
     415            float vf = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
    431416            if (!isfinite(vf)) {
    432417                vf = 1.0;
    … …  
    443428            psFree(readout->variance);
    444429            psFree(readout->covariance);
    445             readout->image  = psMemIncrRefCounter(output->image);
    446             readout->mask   = psMemIncrRefCounter(output->mask);
    447             readout->variance = psMemIncrRefCounter(output->variance);
    448             readout->covariance = psImageCovarianceTruncate(output->covariance, COVAR_FRAC);
    449         } else {
    450             // Fake the convolution
    451             psRegion *region = psRegionAlloc(0, readout->image->numCols - 1, 0, readout->image->numRows - 1);
    452             psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
    453                              PS_DATA_REGION | PS_META_DUPLICATE_OK, "Fake subtraction region", region);
    454             psFree(region);
    455             pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(FAKE_SIZE, 0, penalty,
    456                                                                      PM_SUBTRACTION_MODE_1);
    457             // Set solution to delta function
    458             kernels->solution1 = psVectorAlloc(kernels->num + 2, PS_TYPE_F64);
    459             psVectorInit(kernels->solution1, 0.0);
    460             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    461             kernels->solution1->data.F64[normIndex] = 1.0;
    462             psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
    463                              PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Fake subtraction kernel", kernels);
    464             psFree(kernels);
    465         }
    466 
     430            readout->image  = psMemIncrRefCounter(conv->image);
     431            readout->mask   = psMemIncrRefCounter(conv->mask);
     432            readout->variance = psMemIncrRefCounter(conv->variance);
     433            readout->covariance = psImageCovarianceTruncate(conv->covariance, COVAR_FRAC);
    467434#ifdef TESTING
    468         // Write convolution kernel
     435        }
     436#endif
     437
     438        // Extract the regions and solutions used in the image matching
     439        // This stops them from being freed when we iterate back up the FPA
     440        psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions
    469441        {
    470             const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
    471             assert(outName);
    472 
    473             psString filename = NULL;   // Output filename
    474             psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
    475             psString resolved = pmConfigConvertFilename(filename, config, true, false); // Resolved filename
    476             psFree(filename);
    477             psFits *fits = psFitsOpen(resolved, "w"); // FITS file for subtraction kernel
    478             psFree(resolved);
    479             pmReadoutWriteSubtractionKernels(output, fits);
    480             psFitsClose(fits);
    481         }
    482     }
    483 #endif
    484 
    485     readout->analysis = psMetadataCopy(readout->analysis, output->analysis);
    486 
    487 // Extract the regions and solutions used in the image matching
    488 // This stops them from being freed when we iterate back up the FPA
    489     *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
    490     {
    491         psString regex = NULL;          // Regular expression
    492         psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
    493         psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
    494         psFree(regex);
    495         psMetadataItem *item = NULL;// Item from iteration
    496         while ((item = psMetadataGetAndIncrement(iter))) {
    497             assert(item->type == PS_DATA_REGION);
    498             *regions = psArrayAdd(*regions, ARRAY_BUFFER, item->data.V);
    499         }
    500         psFree(iter);
    501     }
    502     *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels
    503     {
    504         psString regex = NULL;          // Regular expression
    505         psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
    506         psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
    507         psFree(regex);
    508         psMetadataItem *item = NULL;// Item from iteration
    509         while ((item = psMetadataGetAndIncrement(iter))) {
    510             assert(item->type == PS_DATA_UNKNOWN);
    511             // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
    512             // images by scans.
    513             pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
    514             kernel->numCols = readout->image->numCols;
    515             kernel->numRows = readout->image->numRows;
    516 
    517             *kernels = psArrayAdd(*kernels, ARRAY_BUFFER, kernel);
    518         }
    519         psFree(iter);
    520     }
    521     assert((*regions)->n == (*kernels)->n);
    522 
    523     // Record chi^2
    524     {
    525         *chi2 = 0.0;
    526         int num = 0;                    // Number of measurements of chi^2
    527         psString regex = NULL;          // Regular expression
    528         psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
    529         psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
    530         psFree(regex);
    531         psMetadataItem *item = NULL;// Item from iteration
    532         while ((item = psMetadataGetAndIncrement(iter))) {
    533             assert(item->type == PS_DATA_UNKNOWN);
    534             pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
    535             *chi2 += kernels->mean;
    536             num++;
    537         }
    538         psFree(iter);
    539         *chi2 /= psImageCovarianceFactor(readout->covariance) * num;
    540     }
    541 
    542     // Reject image completely if the maximum deconvolution fraction exceeds the limit
    543     float deconv = psMetadataLookupF32(NULL, output->analysis,
    544                                        PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Maximum deconvolution fraction
    545     if (deconv > deconvLimit) {
    546         psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
    547                   deconv, deconvLimit);
    548         psFree(output);
    549         return NULL;
     442            psString regex = NULL;          // Regular expression
     443            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
     444            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
     445            psFree(regex);
     446            psMetadataItem *item = NULL;// Item from iteration
     447            while ((item = psMetadataGetAndIncrement(iter))) {
     448                assert(item->type == PS_DATA_REGION);
     449                regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V);
     450            }
     451            psFree(iter);
     452        }
     453        psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels
     454        {
     455            psString regex = NULL;          // Regular expression
     456            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
     457            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
     458            psFree(regex);
     459            psMetadataItem *item = NULL;// Item from iteration
     460            while ((item = psMetadataGetAndIncrement(iter))) {
     461                assert(item->type == PS_DATA_UNKNOWN);
     462                // Set the normalisation dimensions, since these will be otherwise unavailable when reading
     463                // the images by scans.
     464                pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
     465                kernel->numCols = readout->image->numCols;
     466                kernel->numRows = readout->image->numRows;
     467
     468                kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
     469            }
     470            psFree(iter);
     471        }
     472        psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match");
     473
     474        // Record chi^2
     475        {
     476            double sum = 0.0;           // Sum of chi^2
     477            int num = 0;                // Number of measurements of chi^2
     478            psString regex = NULL;      // Regular expression
     479            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
     480            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
     481            psFree(regex);
     482            psMetadataItem *item = NULL;// Item from iteration
     483            while ((item = psMetadataGetAndIncrement(iter))) {
     484                assert(item->type == PS_DATA_UNKNOWN);
     485                pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
     486                sum += kernels->mean;
     487                num++;
     488            }
     489            psFree(iter);
     490            options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
     491        }
     492
     493        // Reject image completely if the maximum deconvolution fraction exceeds the limit
     494        float deconv = psMetadataLookupF32(NULL, conv->analysis,
     495                                           PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
     496        if (deconv > deconvLimit) {
     497            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
     498                      deconv, deconvLimit);
     499            psFree(conv);
     500            return NULL;
     501        }
     502
     503        readout->analysis = psMetadataCopy(readout->analysis, conv->analysis);
     504
     505        psFree(conv);
     506    } else {
     507        // Match the normalisation
     508        float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
     509        psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
     510        psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
    550511    }
    551512
    552513    // Ensure the background value is zero
    553514    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
     515    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    554516    if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
    555517        psWarning("Can't measure background for image.");
    … …  
    565527    if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
    566528        psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image.");
    567         psFree(output);
     529        psFree(rng);
     530        psFree(bg);
    568531        return false;
    569532    }
    570     *weighting = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
    571                         psImageCovarianceFactor(readout->covariance));
     533    options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
     534                                                  psImageCovarianceFactor(readout->covariance));
    572535    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
    573                      "Weighting by 1/noise^2 for stack", *weighting);
    574 
     536                     "Weighting by 1/noise^2 for stack", options->weightings->data.F32[index]);
     537
     538    psFree(rng);
    575539    psFree(bg);
    576 
    577 #if 0
    578 #define RADIUS 10                       // Radius of photometry
    579 #define MIN_ERR 0.05                    // Minimum photometric error, mag
    580 #define MAX_MAG -13                     // Maximum magnitude for source
    581 
    582     // Ensure the normalisation is correct
    583     // XXX Ideally, would like to do proper PSF-fit photometry, but will settle for aperture photometry
    584     int numSources = sources->n;        // Number of sources
    585     psVector *ratio = psVectorAlloc(numSources, PS_TYPE_F32); // Flux ratios for sources
    586     psVector *ratioMask = psVectorAlloc(numSources, PS_TYPE_MASK); // Mask for flux ratios
    587     psVectorInit(ratioMask, 0xFF);
    588     psImage *image = readout->image;    // Image of interest
    589     psImage *mask = readout->mask;      // Mask of interest
    590     int numCols = image->numCols, numRows = image->numRows; // Size of image
    591     for (int i = 0; i < numSources; i++) {
    592         pmSource *source = sources->data[i]; // Source of interest
    593         if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) ||
    594             source->errMag > MIN_ERR || source->psfMag > MAX_MAG) {
    595             continue;
    596         }
    597 
    598         float xSrc, ySrc;              // Source coordinates
    599         coordsFromSource(&xSrc, &ySrc, source);
    600         int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x
    601         int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y
    602         int numPix = 0;                 // Number of pixels
    603         float sum = 0.0;                // Sum of pixels
    604         for (int y = yMin; y <= yMax; y++) {
    605             for (int x = xMin; x <= xMax; x++) {
    606                 if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) {
    607                     continue;
    608                 }
    609                 sum += image->data.F32[y][x];
    610                 numPix++;
    611             }
    612         }
    613         if (sum >= 0 && numPix > 0) {
    614             float mag = -2.5 * log10(sum * M_PI * PS_SQR(RADIUS) / numPix); // Instrumental magnitude
    615             ratio->data.F32[i] = mag - source->psfMag;
    616             ratioMask->data.PS_TYPE_MASK_DATA[i] = 0;
    617         }
    618     }
    619 
    620     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
    621     if (!psVectorStats(stats, ratio, NULL, ratioMask, 0xFF)) {
    622         psWarning("Unable to measure normalisation --- assuming correct.");
    623     } else {
    624         psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image by %f (+/- %f) mag\n",
    625                  stats->robustMedian, stats->robustStdev);
    626         float norm = powf(10.0, -0.4 * stats->robustMedian); // Normalisation to apply
    627         psBinaryOp(image, image, "*", psScalarAlloc(norm, PS_TYPE_F32));
    628     }
    629     psFree(stats);
    630     psFree(ratio);
    631     psFree(ratioMask);
    632 #endif
    633540
    634541#ifdef TESTING
    … …  
    636543        pmHDU *hdu = pmHDUFromCell(readout->parent);
    637544        psString name = NULL;
    638         psStringAppend(&name, "convolved_%03d.fits", numInput);
    639         pmStackVisualPlotTestImage(output->image, name);
     545        psStringAppend(&name, "convolved_%03d.fits", index);
     546        pmStackVisualPlotTestImage(readout->image, name);
    640547        psFits *fits = psFitsOpen(name, "w");
    641548        psFree(name);
    642         psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
     549        psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
    643550        psFitsClose(fits);
    644551    }
    645552#endif
    646 
    647     psFree(output);
    648553
    649554    return true;

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackOptions.c

    r23352 r23594  
    1818    psFree(options->inputMask);
    1919    psFree(options->sourceLists);
     20    psFree(options->norm);
    2021    psFree(options->cells);
    21     psFree(options->subKernels);
    22     psFree(options->subRegions);
     22    psFree(options->kernels);
     23    psFree(options->regions);
    2324    psFree(options->matchChi2);
    2425    psFree(options->weightings);
    … …  
    3536    psMemSetDeallocator(options, (psFreeFunc)stackOptionsFree);
    3637
     38    options->convolve = true;
    3739    options->stats = NULL;
    3840    options->statsFile = NULL;
    … …  
    4446    options->inputMask = NULL;
    4547    options->sourceLists = NULL;
     48    options->norm = NULL;
    4649    options->cells = NULL;
    47     options->subKernels = NULL;
    48     options->subRegions = NULL;
     50    options->kernels = NULL;
     51    options->regions = NULL;
    4952    options->numCols = 0;
    5053    options->numRows = 0;

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackOptions.h

    r23352 r23594  
    88typedef struct {
    99    // Setup
     10    bool convolve;                      // Convolve images?
    1011    psMetadata *stats;                  // Statistics for output
    1112    FILE *statsFile;                    // File to which to write statistics
    … …  
    1718    psVector *inputMask;                // Mask for inputs
    1819    psArray *sourceLists;               // Individual lists of sources for matching
     20    psVector *norm;                     // Normalisation for each image
    1921    // Convolve
    2022    psArray *cells;                     // Cells for convolved images --- a handle for reading again
    21     psArray *subKernels;                // Subtraction kernels --- required in the stacking
    22     psArray *subRegions;                // Subtraction regions --- required in the stacking
     23    psArray *kernels;                   // PSF-matching kernels --- required in the stacking
     24    psArray *regions;                   // PSF-matching regions --- required in the stacking
    2325    int numCols, numRows;               // Size of image
    2426    psVector *matchChi2;                // chi^2 for stamps from matching

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackPhotometry.c

    r23352 r23594  
    7474        return false;
    7575    }
    76     psFree(photView);
    7776
    7877    if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
    … …  
    9392                         "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
    9493    }
     94    psFree(photView);
     95
    9596    psLogMsg("ppStack", PS_LOG_INFO, "Time to do photometry: %f sec", psTimerClear("PPSTACK_PHOT"));
    9697

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackPrepare.c

    r23352 r23594  
    127127    psVectorInit(options->inputMask, 0);
    128128
    129     if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
    130         pmFPAfileActivate(config->files, false, NULL);
    131         ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
    132         pmFPAview *view = ppStackFilesIterateDown(config);
    133         if (!view) {
    134             return false;
    135         }
    136 
    137         // Generate list of PSFs
    138         psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
    139                                                                "^PPSTACK.INPUT$"); // Iterator
    140         psMetadataItem *fileItem; // Item from iteration
    141         psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
    142         int numCols = 0, numRows = 0;   // Size of image
    143         int index = 0;                  // Index for file
    144         while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
    145             assert(fileItem->type == PS_DATA_UNKNOWN);
    146             pmFPAfile *inputFile = fileItem->data.V; // An input file
     129    pmFPAfileActivate(config->files, false, NULL);
     130    ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
     131    pmFPAview *view = ppStackFilesIterateDown(config);
     132    if (!view) {
     133        return false;
     134    }
     135
     136    psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, "^PPSTACK.INPUT$");
     137    psMetadataItem *fileItem; // Item from iteration
     138    psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
     139    int numCols = 0, numRows = 0;   // Size of image
     140    int index = 0;                  // Index for file
     141    while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
     142        assert(fileItem->type == PS_DATA_UNKNOWN);
     143        pmFPAfile *inputFile = fileItem->data.V; // An input file
     144
     145        // Get list of PSFs, to determine target PSF
     146        if (options->convolve) {
    147147            pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
    148148            pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
    … …  
    173173                numRows = naxis2;
    174174            }
    175 
    176             pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
    177             psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
    178             if (!sources) {
    179                 psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
    180                 return NULL;
    181             }
    182             options->sourceLists->data[index] = psMemIncrRefCounter(sources);
    183 
    184             // Re-do photometry if we don't trust the source lists
    185             if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
    186                 psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
    187                 pmFPAfileActivate(config->files, false, NULL);
    188                 ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index);
    189                 pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
    190                 pmFPAview *photView = ppStackFilesIterateDown(config);
    191                 if (!photView) {
    192                     psFree(view);
    193                     return false;
    194                 }
    195 
    196                 pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
    197 
    198                 if (!ppStackInputPhotometer(ro, sources, config)) {
    199                     psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
    200                     psFree(view);
    201                     psFree(photView);
    202                     return false;
    203                 }
    204 
     175        }
     176
     177
     178        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
     179        psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
     180        if (!sources) {
     181            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
     182            return NULL;
     183        }
     184        options->sourceLists->data[index] = psMemIncrRefCounter(sources);
     185
     186        // Re-do photometry if we don't trust the source lists
     187        if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
     188            psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
     189            pmFPAfileActivate(config->files, false, NULL);
     190            ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index);
     191            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
     192            pmFPAview *photView = ppStackFilesIterateDown(config);
     193            if (!photView) {
     194                psFree(view);
     195                return false;
     196            }
     197
     198            pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
     199
     200            if (!ppStackInputPhotometer(ro, sources, config)) {
     201                psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
     202                psFree(view);
    205203                psFree(photView);
    206                 if (!ppStackFilesIterateUp(config)) {
    207                     psFree(view);
    208                     return false;
    209                 }
    210                 pmFPAfileActivate(config->files, false, NULL);
    211                 ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
    212             }
    213 
    214             index++;
    215         }
    216         psFree(fileIter);
    217 
    218         // Zero point calibration
    219         options->sumExposure = ppStackSourcesTransparency(options->sourceLists, options->inputMask,
    220                                                           view, config);
    221         if (!isfinite(options->sumExposure) || options->sumExposure <= 0) {
    222             psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
    223             psFree(view);
    224             return false;
    225         }
    226 
    227         // Generate target PSF
     204                return false;
     205            }
     206
     207            psFree(photView);
     208            if (!ppStackFilesIterateUp(config)) {
     209                psFree(view);
     210                return false;
     211            }
     212            pmFPAfileActivate(config->files, false, NULL);
     213            ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
     214        }
     215
     216        index++;
     217    }
     218    psFree(fileIter);
     219
     220    // Generate target PSF
     221    if (options->convolve) {
    228222        options->psf = ppStackPSF(config, numCols, numRows, psfs, options->inputMask);
    229223        psFree(psfs);
    … …  
    240234                         "Target PSF", options->psf);
    241235        outChip->data_exists = true;
    242 
     236    }
     237
     238    // Zero point calibration
     239    if (!ppStackSourcesTransparency(options, view, config)) {
     240        psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
    243241        psFree(view);
    244 
    245         if (!ppStackFilesIterateUp(config)) {
    246             return false;
    247         }
     242        return false;
     243    }
     244    psFree(view);
     245
     246    if (!ppStackFilesIterateUp(config)) {
     247        return false;
    248248    }
    249249

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackReadout.c

    r23352 r23594  
    150150
    151151    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
    152                         combineRej, combineSys, combineDiscard, true, useVariance, safe)) {
     152                        combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) {
    153153        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
    154154        psFree(stack);
    … …  
    185185    assert(outRO);
    186186    assert(readouts);
    187     assert(rejected);
    188     assert(readouts->n == rejected->n);
     187    assert(!rejected || readouts->n == rejected->n);
    189188    assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK);
    190189    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
    … …  
    197196
    198197    bool mdok;                          // Status of MD lookup
     198    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
     199    float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
     200    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
     201    float combineDiscard = psMetadataLookupF32(NULL, recipe, "COMBINE.DISCARD"); // Olympic discard fraction
    199202    bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
     203    bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
    200204
    201205    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
    … …  
    207211    psArray *stack = psArrayAlloc(num); // Array for stacking
    208212
    209     int numGood = num;                  // Number of good inputs: images that haven't been completely rejected
     213    bool entire = true;                 // Combine entire image?
     214    if (rejected) {
     215        // We have rejection from a previous combination: combine without flagging pixels to inspect
     216        entire = false;
     217        safe = false;
     218        iter = 0;
     219        combineRej = NAN;
     220        combineSys = NAN;
     221    }
     222
    210223    for (int i = 0; i < num; i++) {
    211224        pmReadout *ro = readouts->data[i];
    212         if (!ro || !rejected->data[i] || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    213             // Image completely rejected
    214             numGood--;
     225        if (mask->data.U8[i] & (PPSTACK_MASK_REJECT | PPSTACK_MASK_BAD)) {
     226            // Image completely rejected since previous combination
     227            entire = true;
     228            continue;
     229        } else if (mask->data.U8[i]) {
     230            // Image completely rejected before original combination
    215231            continue;
    216232        }
    … …  
    234250        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i],
    235251                                             addVariance ? addVariance->data.F32[i] : NAN);
    236         data->reject = psMemIncrRefCounter(rejected->data[i]);
     252        data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL;
    237253        stack->data[i] = data;
    238254    }
    239255
    240     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, 0, NAN, NAN, NAN,
    241                         numGood != num, useVariance, false)) {
     256    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0,
     257                        iter, combineRej, combineSys, combineDiscard,
     258                        entire, useVariance, safe, !rejected)) {
    242259        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
    243260        psFree(stack);

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackReject.c

    r23352 r23594  
    1414    psAssert(options, "Require options");
    1515    psAssert(config, "Require configuration");
     16
     17    if (!options->convolve) {
     18        // No need to do complicated rejection when we haven't convolved
     19        return true;
     20    }
    1621
    1722    int num = options->num;             // Number of inputs
    … …  
    8691
    8792        psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows,
    88                                          threshold, poorFrac, stride, options->subRegions->data[i],
    89                                          options->subKernels->data[i]); // Rejected pixels
     93                                         threshold, poorFrac, stride, options->regions->data[i],
     94                                         options->kernels->data[i]); // Rejected pixels
    9095
    9196#ifdef TESTING
    … …  
    141146
    142147    psFree(options->inspect); options->inspect = NULL;
    143     psFree(options->subKernels); options->subKernels = NULL;
    144     psFree(options->subRegions); options->subRegions = NULL;
     148    psFree(options->kernels); options->kernels = NULL;
     149    psFree(options->regions); options->regions = NULL;
    145150
    146151    if (numRejected >= num - 1) {

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackSetup.c

    r23352 r23594  
    1111#include "ppStackLoop.h"
    1212
     13#define BUFFER 16                       // Buffer for name array
     14
     15// Generate an array of input filenames
     16static psArray *stackNameArray(const ppStackOptions *options, // Stack options
     17                               pmConfig *config, // Configuration
     18                               const char *name // Name of file
     19    )
     20{
     21    psAssert(config, "Require configuration");
     22    psAssert(config, "Require file name");
     23
     24    psArray *array = psArrayAllocEmpty(options->num); // Array with filenames
     25    pmFPAview *view = pmFPAviewAlloc(0);// View to readout
     26    view->chip = view->cell = view->readout = -1;
     27
     28    for (int i = 0; i < options->num; i++) {
     29        pmFPAfile *file = pmFPAfileSelectSingle(config->files, name, i); // File of interest
     30
     31        psString filename = pmFPAfileName(file, view, config); // Filename of interest
     32        psAssert(filename, "Can't determine filename");
     33        psArrayAdd(array, array->n, filename);
     34        psFree(filename);               // Drop reference
     35    }
     36
     37    return array;
     38}
     39
     40
    1341bool ppStackSetup(ppStackOptions *options, pmConfig *config)
    1442{
    … …  
    1846    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    1947    psAssert(recipe, "We've thrown an error on this before.");
     48
     49    options->convolve = psMetadataLookupBool(NULL, config->arguments, "CONVOLVE"); // Convolve images?
     50    if (!psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
     51        psWarning("No PSFs provided --- unable to convolve to common PSF.");
     52        options->convolve = false;
     53    }
    2054
    2155    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
    … …  
    3670    }
    3771
    38     const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images
    39     if (!tempDir) {
    40         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe");
    41         return false;
     72    // Generate temporary names for convolved images
     73    if (options->convolve) {
     74        const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images
     75        if (!tempDir) {
     76            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe");
     77            return false;
     78        }
     79
     80        psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments,
     81                                                               "OUTPUT")); // Name for temporary files
     82        const char *tempName = basename(outputName);
     83        if (!tempName) {
     84            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files.");
     85            psFree(outputName);
     86            return false;
     87        }
     88
     89        const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for images
     90        const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for masks
     91        const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for var maps
     92        if (!tempImage || !tempMask || !tempVariance) {
     93            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     94                    "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
     95            psFree(outputName);
     96            return false;
     97        }
     98
     99        options->imageNames = psArrayAlloc(num);
     100        options->maskNames = psArrayAlloc(num);
     101        options->varianceNames = psArrayAlloc(num);
     102        for (int i = 0; i < num; i++) {
     103            psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
     104            psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage);
     105            psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask);
     106            psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance);
     107            psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName);
     108            options->imageNames->data[i] = imageName;
     109            options->maskNames->data[i] = maskName;
     110            options->varianceNames->data[i] = varianceName;
     111        }
     112        psFree(outputName);
     113    } else {
     114        options->imageNames = stackNameArray(options, config, "PPSTACK.INPUT");
     115        options->maskNames = stackNameArray(options, config, "PPSTACK.INPUT.MASK");
     116        options->varianceNames = stackNameArray(options, config, "PPSTACK.INPUT.VARIANCE");
    42117    }
    43 
    44     psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments,
    45                                                            "OUTPUT")); // Name for temporary files
    46     const char *tempName = basename(outputName);
    47     if (!tempName) {
    48         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files.");
    49         psFree(outputName);
    50         return false;
    51     }
    52 
    53     const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for temporary images
    54     const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for temporary masks
    55     const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for temp var maps
    56     if (!tempImage || !tempMask || !tempVariance) {
    57         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    58                 "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
    59         psFree(outputName);
    60         return false;
    61     }
    62 
    63     options->imageNames = psArrayAlloc(num);
    64     options->maskNames = psArrayAlloc(num);
    65     options->varianceNames = psArrayAlloc(num);
    66     for (int i = 0; i < num; i++) {
    67         psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
    68         psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage);
    69         psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask);
    70         psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance);
    71         psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName);
    72         options->imageNames->data[i] = imageName;
    73         options->maskNames->data[i] = maskName;
    74         options->varianceNames->data[i] = varianceName;
    75     }
    76     psFree(outputName);
    77118
    78119    if (!pmConfigMaskSetBits(NULL, NULL, config)) {

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackSources.c

    r23352 r23594  
    99//#define TESTING                         // Enable debugging output
    1010
     11#ifdef TESTING
    1112// Size of fake image; set by hand because it's trouble to get it from other places
    1213#define FAKE_COLS 4861
    1314#define FAKE_ROWS 4913
     15#endif
    1416
    1517#ifdef TESTING
    … …  
    5355
    5456
    55 float ppStackSourcesTransparency(const psArray *sourceLists, psVector *inputMask,
    56                                  const pmFPAview *view, const pmConfig *config)
     57bool ppStackSourcesTransparency(ppStackOptions *options, const pmFPAview *view, const pmConfig *config)
    5758{
    58     PS_ASSERT_ARRAY_NON_NULL(sourceLists, NAN);
    59     PS_ASSERT_VECTOR_NON_NULL(inputMask, NAN);
    60     PS_ASSERT_VECTOR_TYPE(inputMask, PS_TYPE_U8, NAN);
    61     PS_ASSERT_VECTOR_SIZE(inputMask, sourceLists->n, NAN);
    62     PS_ASSERT_PTR_NON_NULL(view, NAN);
    63     PS_ASSERT_PTR_NON_NULL(config, NAN);
     59    PS_ASSERT_PTR_NON_NULL(options, false);
     60    PS_ASSERT_PTR_NON_NULL(view, false);
     61    PS_ASSERT_PTR_NON_NULL(config, false);
     62
     63    psArray *sourceLists = options->sourceLists; // Source lists for each input
     64    psVector *inputMask = options->inputMask; // Mask for inputs
     65
     66    PS_ASSERT_ARRAY_NON_NULL(sourceLists, false);
     67    PS_ASSERT_VECTOR_NON_NULL(inputMask, false);
     68    PS_ASSERT_VECTOR_TYPE(inputMask, PS_TYPE_U8, false);
     69    PS_ASSERT_VECTOR_SIZE(inputMask, sourceLists->n, false);
    6470
    6571#if defined(TESTING) && 0
    … …  
    100106    if (!airmassZP) {
    101107        psError(PS_ERR_UNKNOWN, false, "Unable to find ZP.AIRMASS in recipe.");
    102         return NAN;
     108        return false;
    103109    }
    104110
    … …  
    139145                    exptime, airmass, expFilter);
    140146            psFree(zp);
    141             return NAN;
     147            return false;
    142148        }
    143149
    … …  
    149155                        "Unable to find airmass term (ZP.AIRMASS) for filter %s", filter);
    150156                psFree(zp);
    151                 return NAN;
     157                return false;
    152158            }
    153159        } else if (strcmp(filter, expFilter) != 0) {
    154160            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Filters don't match: %s vs %s", filter, expFilter);
    155161            psFree(zp);
    156             return NAN;
     162            return false;
    157163        }
    158164
    … …  
    160166        sumExpTime += exptime;
    161167    }
     168
     169    options->sumExposure = sumExpTime;
    162170
    163171    psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
    … …  
    165173        psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
    166174        psFree(zp);
    167         return NAN;
     175        return false;
    168176    }
    169177
    … …  
    177185    if (!trans) {
    178186        psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
    179         return NAN;
     187        return false;
    180188    }
    181189
    … …  
    186194    for (int i = 0; i < trans->n; i++) {
    187195        if (!isfinite(trans->data.F32[i])) {
    188             inputMask->data.U8[i] = PPSTACK_MASK_CAL;
     196            inputMask->data.U8[i] |= PPSTACK_MASK_CAL;
    189197        }
    190198    }
    … …  
    223231    // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
    224232    // We don't need to know the magnitude zero point for the filter, since it cancels out
     233    options->norm = psVectorAlloc(num, PS_TYPE_F32);
    225234    for (int i = 0; i < num; i++) {
    226235        if (!isfinite(trans->data.F32[i])) {
    … …  
    229238        psArray *sources = sourceLists->data[i]; // Sources of interest
    230239        float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];
     240        options->norm->data.F32[i] = magCorr;
    231241        psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr);
    232242
    … …  
    248258            psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
    249259            psFree(zp);
    250             return NAN;
     260            return false;
    251261        }
    252262        psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
    … …  
    259269#endif
    260270
    261     return sumExpTime;
     271    return true;
    262272}
Note: See TracChangeset for help on using the changeset viewer.