IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 17, 2009, 12:08:50 PM (17 years ago)
Author:
beaumont
Message:

merged with head

Location:
branches/cnb_branches/cnb_branch_20090301
Files:
14 edited
13 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

    r19337 r23352  
    11bin_PROGRAMS = ppStack
    22
    3 ppStack_CFLAGS  = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) $(PPSTATS_CFLAGS) $(PPSTACK_CFLAGS)
     3# PPSTACK_VERSION=`if [ -e ../../VERSION ]; then cat ../../VERSION; else svnversion; fi`
     4# PPSTACK_BRANCH=`if [ -e ../../BRANCH ]; then cat ../../BRANCH; else svn info | sed -n -e '/URL:/ h' -e '/Repository Root:/ { x; H; x; s|Repository Root: \(.*\)\nURL: \1\(.*\)|\2| ; s|^/|| ; s|/[a-zA-Z]*/src.*|| ; p }'; fi`
     5# PPSTACK_SOURCE=`if [ -e ../../SOURCE ]; then cat ../../SOURCE; else svn info | sed -n -e 's/Repository UUID: // p'; fi`
     6#
     7# # Force recompilation of ppStackVersion.c, since it gets the version information
     8# ppStackVersion.c: FORCE
     9#       touch ppStackVersion.c
     10# FORCE: ;
     11
     12ppStack_CFLAGS  = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) $(PPSTATS_CFLAGS) $(PPSTACK_CFLAGS) -DPPSTACK_VERSION=$(SVN_VERSION) -DPPSTACK_BRANCH=$(SVN_BRANCH) -DPPSTACK_SOURCE=$(SVN_SOURCE)
    413ppStack_LDFLAGS = $(PSLIB_LIBS)   $(PSMODULE_LIBS)   $(PSPHOT_LIBS)   $(PPSTATS_LIBS)   $(PPSTACK_LIBS)
    514
    … …  
    817        ppStackArguments.c      \
    918        ppStackCamera.c         \
     19        ppStackFiles.c          \
    1020        ppStackLoop.c           \
    1121        ppStackPSF.c            \
    1222        ppStackReadout.c        \
    13         ppStackPhotometry.c     \
    1423        ppStackVersion.c        \
    1524        ppStackMatch.c          \
    1625        ppStackSources.c        \
    17         ppStackThread.c
     26        ppStackThread.c         \
     27        ppStackOptions.c        \
     28        ppStackSetup.c          \
     29        ppStackPrepare.c        \
     30        ppStackConvolve.c       \
     31        ppStackCombineInitial.c \
     32        ppStackReject.c         \
     33        ppStackCombineFinal.c   \
     34        ppStackCleanup.c        \
     35        ppStackPhotometry.c     \
     36        ppStackFinish.c
    1837
    1938noinst_HEADERS =                \
    20         ppStack.h
     39        ppStack.h               \
     40        ppStackLoop.h           \
     41        ppStackOptions.h        \
     42        ppStackThread.h
    2143
    2244CLEANFILES = *~

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

    r21258 r23352  
    99
    1010#include "ppStack.h"
     11#include "ppStackLoop.h"
    1112
    1213#define TIMER_NAME "PPSTACK"            // Name of timer
    … …  
    3334        goto die;
    3435    }
     36
     37    ppStackVersionPrint();
    3538
    3639    if (!pmModelClassInit()) {
    … …  
    7477    pmConfigDone();
    7578    psLibFinalize();
    76 
     79    pmVisualClose();
    7780    exit(exitValue);
    7881}

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

    r21477 r23352  
    1 #ifndef PP_STACK_H
    2 #define PP_STACK_H
     1#ifndef PPSTACK_H
     2#define PPSTACK_H
    33
    44#define PPSTACK_RECIPE "PPSTACK"        // Name of the recipe
    … …  
    1010// Mask values for inputs
    1111typedef enum {
    12     PPSTACK_MASK_MATCH  = 0x01,         // PSF-matching failed
    13     PPSTACK_MASK_CHI2   = 0x02,         // Chi^2 too deviant
    14     PPSTACK_MASK_REJECT = 0x04,         // Rejection failed
    15     PPSTACK_MASK_BAD    = 0x08,         // Bad image (too many pixels rejected)
     12    PPSTACK_MASK_CAL    = 0x01,         // Photometric calibration failed
     13    PPSTACK_MASK_MATCH  = 0x02,         // PSF-matching failed
     14    PPSTACK_MASK_CHI2   = 0x04,         // Chi^2 too deviant
     15    PPSTACK_MASK_REJECT = 0x08,         // Rejection failed
     16    PPSTACK_MASK_BAD    = 0x10,         // Bad image (too many pixels rejected)
    1617    PPSTACK_MASK_ALL    = 0xff          // All errors
    1718} ppStackMask;
    1819
    19 // Thread for stacking chunks
    20 //
    21 // Each input file contributes a readout, into which is read a chunk from that file
    22 typedef struct {
    23     psArray *readouts;                  // Input readouts to read and stack
    24     bool read;                          // Has the scan been read?
    25     bool busy;                          // Is the scan being processed?
    26     int firstScan;                      // First row of the chunk to be read for this group
    27     int lastScan;                       // Last row of the chunk to be read for this group
    28 } ppStackThread;
     20// List of files
     21typedef enum {
     22    PPSTACK_FILES_PREPARE,              // Files for preparation
     23    PPSTACK_FILES_CONVOLVE,             // Files for convolution
     24    PPSTACK_FILES_COMBINE,              // Files for combination
     25    PPSTACK_FILES_PHOT                  // Files for photometry
     26} ppStackFileList;
    2927
    30 // Allocator
    31 ppStackThread *ppStackThreadAlloc(psArray *readouts // Inputs readouts to read and stack
    32     );
    33 
    34 // Data for threads
    35 typedef struct {
    36     psArray *threads;                   // Threads doing stacking
    37     int lastScan;                       // Last row that's been read
    38     psArray *imageFits;                 // FITS file pointers for images
    39     psArray *maskFits;                  // FITS file pointers for masks
    40     psArray *varianceFits;                // FITS file pointers for variances
    41 } ppStackThreadData;
    42 
    43 // Set up thread data
    44 ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, // Array of input cells
    45                                           const psArray *imageNames, // Names of images to read
    46                                           const psArray *maskNames, // Names of masks to read
    47                                           const psArray *varianceNames, // Names of variance maps to read
    48                                           const psArray *covariances, // Covariance matrices (already read)
    49                                           const pmConfig *config // Configuration
    50     );
    51 
    52 // Read chunk into the first available file thread
    53 ppStackThread *ppStackThreadRead(bool *status, // Status of read
    54                                  ppStackThreadData *stack, // Stacks available for reading
    55                                  pmConfig *config, // Configuration
    56                                  int numChunk, // Chunk number (only for interest)
    57                                  int overlap // Overlap between subsequent scans
    58     );
    59 
    60 // Initialise the threads
    61 void ppStackThreadInit(void);
    6228
    6329// Setup command-line arguments
    … …  
    7440    );
    7541
    76 // Loop over the inputs, doing the combination
    77 bool ppStackLoop(pmConfig *config       // Configuration
    78     );
    79 
    8042// Determine target PSF for input images
    8143pmPSF *ppStackPSF(const pmConfig *config, // Configuration
    8244                  int numCols, int numRows, // Size of image
    83                   const psArray *psfs   // List of input PSFs
     45                  const psArray *psfs,  // List of input PSFs
     46                  const psVector *inputMask // Mask for inputs
    8447    );
    8548
    … …  
    12588    );
    12689
    127 // Perform photometry on stack
    128 bool ppStackPhotometry(pmConfig *config, // Configuration
    129                        const pmReadout *readout, // Readout to be photometered
    130                        const pmFPAview *view // View to readout
    131     );
    132 
    13390// Return software version
    13491psString ppStackVersion(void);
     92
     93/// Return software source
     94psString ppStackSource(void);
    13595
    13696// Return long description of software version
    13797psString ppStackVersionLong(void);
    13898
    139 // Supplement metadata with software version
    140 void ppStackVersionMetadata(psMetadata *metadata // Metadata to supplement
     99// Supplement header with software version
     100bool ppStackVersionHeader(psMetadata *header // Header to supplement
    141101    );
     102
     103/// Print version information
     104void ppStackVersionPrint(void);
    142105
    143106/// Convolve image to match specified seeing
    … …  
    158121/// Corrects the source PSF photometry to a common system.  Return the sum of the exposure times.
    159122float ppStackSourcesTransparency(const psArray *sourceLists, // Sources for each input
     123                                 psVector *inputMask, // Indicates bad input
    160124                                 const pmFPAview *view, // View to readout
    161125                                 const pmConfig *config // Configuration
    162126    );
    163127
     128/// Dump memory debugging information
     129void ppStackMemDump(
     130    const char *name                    ///< Stage name, for inclusion in the output file name
     131    );
     132
     133/// Activate/deactivate a list of files
     134void ppStackFileActivation(
     135    pmConfig *config,                   // Configuration
     136    ppStackFileList list,               // Files to turn on/off
     137    bool state                          // Activation state
     138    );
     139
     140// Activate/deactivate a single element for a list
     141void ppStackFileActivationSingle(
     142    pmConfig *config,                   // Configuration
     143    ppStackFileList list,               // Files to turn on/off
     144    bool state,                         // Activation state
     145    int num                             // Number of file in sequence
     146    );
     147
     148/// Iterate down the hierarchy, loading files
     149///
     150/// We can get away with this simplistic treatment of the FPA hierarchy because we're working on skycells.
     151pmFPAview *ppStackFilesIterateDown(
     152    pmConfig *config                    // Configuration
     153    );
     154
     155/// Iterate up the hierarchy, writing files
     156///
     157/// We can get away with this simplistic treatment of the FPA hierarchy because we're working on skycells.
     158bool ppStackFilesIterateUp(
     159    pmConfig *config                    // Configuration
     160    );
     161
     162/// Write an image to a FITS file
     163bool ppStackWriteImage(
     164    const char *name,                   // Name of image
     165    psMetadata *header,                 // Header
     166    const psImage *image,               // Image
     167    pmConfig *config                    // Configuration
     168    );
     169
     170
    164171#endif

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

    r21525 r23352  
    168168                      "Play safe with small numbers of pixels to combine?", false);
    169169    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-radius", 0, "Radius (pixels) for matching sources", NAN);
    170     psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-iter", 0, "Maximum iterations for zero point", 0);
     170    psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-iter-1", 0, "Maximum iterations for zero point; pass 1", 0);
     171    psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-iter-2", 0, "Maximum iterations for zero point; pass 2", 0);
    171172    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-tol", 0, "Tolerance for zero point iterations", NAN);
    172173    psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-trans-iter", 0, "Iterations for transparency determination", 0);
    173174    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-trans-rej", 0, "Rejection threshold for transparency determination", NAN);
    174175    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-trans-thresh", 0, "Threshold for transparency determination", NAN);
    175     psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-rej", 0, "Rejection threshold for stars", NAN);
     176    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-rej-1", 0, "Rejection threshold for stars; pass 1", NAN);
     177    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-rej-2", 0, "Rejection threshold for stars; pass 2", NAN);
    176178    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-limit", 0, "Limit on star rejection fraction for successful iteration", NAN);
    177     psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys", 0, "Estimated systematic error", NAN);
     179    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys-1", 0, "Estimated systematic error; pass 1", NAN);
     180    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys-2", 0, "Estimated systematic error; pass 2", NAN);
    178181    psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-image", 0, "Suffix for temporary images", NULL);
    179182    psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-mask", 0, "Suffix for temporary masks", NULL);
    … …  
    182185                      "Delete temporary files on completion?", false);
    183186    psMetadataAddS32(arguments, PS_LIST_TAIL, "-threads", 0, "Number of threads to use", 0);
    184     psMetadataAddBool(arguments, PS_LIST_TAIL, "-psphot-visual", 0, "psphot visualisation", 0);
     187    psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "visualisation", false);
    185188
    186189    if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) {
    … …  
    242245
    243246    VALUE_ARG_RECIPE_FLOAT("-zp-radius",       "ZP.RADIUS",       F32);
    244     VALUE_ARG_RECIPE_INT("-zp-iter",           "ZP.ITER",         S32, 0);
     247    VALUE_ARG_RECIPE_INT("-zp-iter-1",         "ZP.ITER.1",       S32, 0);
     248    VALUE_ARG_RECIPE_INT("-zp-iter-2",         "ZP.ITER.2",       S32, 0);
    245249    VALUE_ARG_RECIPE_FLOAT("-zp-tol",          "ZP.TOL",          F32);
    246250    VALUE_ARG_RECIPE_INT("-zp-trans-iter",     "ZP.TRANS.ITER",   S32, 0);
    247251    VALUE_ARG_RECIPE_FLOAT("-zp-trans-rej",    "ZP.TRANS.REJ",    F32);
    248252    VALUE_ARG_RECIPE_FLOAT("-zp-trans-thresh", "ZP.TRANS.THRESH", F32);
    249     VALUE_ARG_RECIPE_FLOAT("-zp-star-rej",     "ZP.STAR.REJ",     F32);
     253    VALUE_ARG_RECIPE_FLOAT("-zp-star-rej-1",   "ZP.STAR.REJ.1",   F32);
     254    VALUE_ARG_RECIPE_FLOAT("-zp-star-rej-2",   "ZP.STAR.REJ.2",   F32);
    250255    VALUE_ARG_RECIPE_FLOAT("-zp-star-limit",   "ZP.STAR.LIMIT",   F32);
    251     VALUE_ARG_RECIPE_FLOAT("-zp-star-sys",     "ZP.STAR.SYS",     F32);
     256    VALUE_ARG_RECIPE_FLOAT("-zp-star-sys-1",   "ZP.STAR.SYS.1",   F32);
     257    VALUE_ARG_RECIPE_FLOAT("-zp-star-sys-2",   "ZP.STAR.SYS.2",   F32);
    252258
    253259    VALUE_ARG_RECIPE_INT("-psf-instances", "PSF.INSTANCES", S32, 0);
    … …  
    255261    VALUE_ARG_RECIPE_INT("-psf-order",     "PSF.ORDER",     S32, 0);
    256262    valueArgRecipeStr(arguments, recipe, "-psf-model", "PSF.MODEL", recipe);
     263
     264    if (psMetadataLookupBool(NULL, arguments, "-visual")) {
     265        pmVisualSetVisual(true);
     266    }
    257267
    258268    if (psMetadataLookupBool(NULL, arguments, "-photometry") ||
    … …  
    290300    if (dump_file) {
    291301        pmConfigCamerasCull(config, NULL);
    292         pmConfigRecipesCull(config, "PPSTACK,PPSUB,PPSTATS,PSPHOT,MASKS");
     302        pmConfigRecipesCull(config, "PPSTACK,PPSUB,PPSTATS,PSPHOT,MASKS,JPEG");
    293303
    294304        pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Input file

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

    r21526 r23352  
    44
    55#include <stdio.h>
    6 #include <unistd.h>
    7 #include <string.h>
    8 #include <libgen.h>
    96#include <pslib.h>
    107#include <psmodules.h>
    11 #include <ppStats.h>
    128
    139#include "ppStack.h"
    14 
    15 //#define TESTING
    16 
    17 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
    18 #define PIXELS_BUFFER 1024              // Initial size of pixel lists
    19 
    20 // Here follows lists of files for activation/deactivation at various stages.  Each must be NULL-terminated.
    21 
    22 // Files required in preparation for convolution
    23 static char *prepareFiles[] = { "PPSTACK.INPUT.PSF", "PPSTACK.INPUT.SOURCES", "PPSTACK.TARGET.PSF", NULL };
    24 
    25 // Files required for the convolution
    26 static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL };
    27 
    28 // Output files for the combination
    29 static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE",
    30                                 "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL };
    31 
    32 // Files for photometry
    33 static char *photFiles[] = { "PSPHOT.INPUT", "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL",
    34                              "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
    35                              "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID",
    36                              "PSPHOT.INPUT.CMF", NULL };
    37 
    38 static void memDump(const char *name)
    39 {
    40     return;
    41 
    42     static int num = 0;                 // Counter, to make files unique and give an idea of sequence
    43 
    44     psString filename = NULL;           // Name of file
    45     psStringAppend(&filename, "memdump_%s_%03d.txt", name, num);
    46     FILE *memFile = fopen(filename, "w");
    47     psFree(filename);
    48 
    49     psMemBlock **leaks = NULL;
    50     int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true);
    51     fprintf(memFile, "# MemBlock Size Source\n");
    52     unsigned long total = 0;            // Total memory used
    53     for (int i = 0; i < numLeaks; i++) {
    54         psMemBlock *mb = leaks[i];
    55         fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize,
    56                 mb->file, mb->lineno);
    57         total += mb->userMemorySize;
    58     }
    59     fclose(memFile);
    60     psFree(leaks);
    61 
    62     fprintf(stderr, "Memdump %s %d: Memory use: %ld, sbrk: %p\n", name, num, total, sbrk(0));
    63     num++;
    64 }
    65 
    66 
    67 
    68 // Activate/deactivate a list of files
    69 static void fileActivation(pmConfig *config, // Configuration
    70                            char **files, // Files to turn on/off
    71                            bool state   // Activation state
    72     )
    73 {
    74     assert(config);
    75     assert(files);
    76 
    77     for (int i = 0; files[i] != NULL; i++) {
    78         pmFPAfileActivate(config->files, state, files[i]);
    79     }
    80     return;
    81 }
    82 
    83 // Activate/deactivate a single element for a list
    84 static void fileActivationSingle(pmConfig *config, // Configuration
    85                                  char **files, // Files to turn on/off
    86                                  bool state,   // Activation state
    87                                  int num // Number of file in sequence
    88                                  )
    89 {
    90     assert(config);
    91     assert(files);
    92 
    93     for (int i = 0; files[i] != NULL; i++) {
    94         pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file
    95     }
    96     return;
    97 }
    98 
    99 // Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells
    100 static pmFPAview *filesIterateDown(pmConfig *config // Configuration
    101                                   )
    102 {
    103     assert(config);
    104 
    105     pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
    106     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    107         return NULL;
    108     }
    109     view->chip = 0;
    110     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    111         return NULL;
    112     }
    113     view->cell = 0;
    114     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    115         return NULL;
    116     }
    117     view->readout = 0;
    118     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    119         return NULL;
    120     }
    121     return view;
    122 }
    123 
    124 // Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells
    125 static bool filesIterateUp(pmConfig *config // Configuration
    126                            )
    127 {
    128     assert(config);
    129 
    130     pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
    131     view->chip = view->cell = view->readout = 0;
    132     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    133         return false;
    134     }
    135     view->readout = -1;
    136     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    137         return false;
    138     }
    139     view->cell = -1;
    140     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    141         return false;
    142     }
    143     view->chip = -1;
    144     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    145         return false;
    146     }
    147     psFree(view);
    148     return true;
    149 }
    150 
    151 // Write an image to a FITS file
    152 static bool writeImage(const char *name, // Name of image
    153                        psMetadata *header, // Header
    154                        const psImage *image, // Image
    155                        pmConfig *config // Configuration
    156                        )
    157 {
    158     assert(name);
    159     assert(image);
    160 
    161     psString resolved = pmConfigConvertFilename(name, config, true, true); // Resolved file name
    162     psFits *fits = psFitsOpen(resolved, "w");
    163     if (!fits) {
    164         psError(PS_ERR_IO, false, "Unable to open FITS file %s to write image.", resolved);
    165         psFree(resolved);
    166         return false;
    167     }
    168     if (!psFitsWriteImage(fits, header, image, 0, NULL)) {
    169         psError(PS_ERR_IO, false, "Unable to write FITS image %s.", resolved);
    170         psFitsClose(fits);
    171         psFree(resolved);
    172         return false;
    173     }
    174     psFitsClose(fits);
    175     psFree(resolved);
    176     return true;
    177 }
    178 
     10#include "ppStackLoop.h"
    17911
    18012bool ppStackLoop(pmConfig *config)
    … …  
    18315
    18416    psTimerStart("PPSTACK_TOTAL");
     17    psTimerStart("PPSTACK_STEPS");
    18518
    186     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    187     psAssert(recipe, "We've thrown an error on this before.");
     19    ppStackOptions *options = ppStackOptionsAlloc(); // Options for stacking
    18820
    189     bool mdok;                          // Status of MD lookup
    190     bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files?
    191     const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images
    192     if (!tempDir) {
    193         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe");
     21    // Setup
     22    psTrace("ppStack", 1, "Setup....\n");
     23    if (!ppStackSetup(options, config)) {
     24        psError(PS_ERR_UNKNOWN, false, "Unable to setup.");
     25        psFree(options);
     26        return false;
     27    }
     28    psLogMsg("ppStack", PS_LOG_INFO, "Stage 0: Setup: %f sec", psTimerClear("PPSTACK_STEPS"));
     29    ppStackMemDump("setup");
     30
     31
     32    // Preparation for stacking
     33    psTrace("ppStack", 1, "Preparation for stacking: merging sources, determining target PSF....\n");
     34    if (!ppStackPrepare(options, config)) {
     35        psError(PS_ERR_UNKNOWN, false, "Unable to prepare for stacking.");
     36        psFree(options);
     37        return false;
     38    }
     39    psLogMsg("ppStack", PS_LOG_INFO, "Stage 1: Load Sources and Generate Target PSF: %f sec",
     40             psTimerClear("PPSTACK_STEPS"));
     41    ppStackMemDump("prepare");
     42
     43
     44    // Convolve inputs
     45    psTrace("ppStack", 1, "Convolving inputs to target PSF....\n");
     46    if (!ppStackConvolve(options, config)) {
     47        psError(PS_ERR_UNKNOWN, false, "Unable to convolve images.");
     48        psFree(options);
     49        return false;
     50    }
     51    psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec",
     52             psTimerClear("PPSTACK_STEPS"));
     53    ppStackMemDump("convolve");
     54
     55
     56    // Start threading
     57    ppStackThreadInit();
     58    ppStackThreadData *stack = ppStackThreadDataSetup(options, config);
     59    if (!stack) {
     60        psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
     61        psFree(options);
    19462        return false;
    19563    }
    19664
    197     psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments, "OUTPUT")); // Name for temporary files
    198     const char *tempName = basename(outputName);
    199     if (!tempName) {
    200         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files.");
     65    // Initial combination
     66    psTrace("ppStack", 1, "Initial stack of convolved images....\n");
     67    if (!ppStackCombineInitial(stack, options, config)) {
     68        psError(PS_ERR_UNKNOWN, false, "Unable to perform initial combination.");
     69        psFree(stack);
     70        psFree(options);
    20171        return false;
    20272    }
     73    psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec",
     74             psTimerClear("PPSTACK_STEPS"));
     75    ppStackMemDump("convolve");
     76    psLogMsg("ppStack", PS_LOG_INFO, "Stage 4: Make Initial Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
    20377
    204     const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for temporary images
    205     const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for temporary masks
    206     const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for temp variance maps
    207     if (!tempImage || !tempMask || !tempVariance) {
    208         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    209                 "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
     78
     79    // Pixel rejection
     80    psTrace("ppStack", 1, "Reject pixels....\n");
     81    if (!ppStackReject(options, config)) {
     82        psError(PS_ERR_UNKNOWN, false, "Unable to reject pixels.");
     83        psFree(stack);
     84        psFree(options);
    21085        return false;
    21186    }
     87    psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Pixel Rejection: %f sec", psTimerClear("PPSTACK_STEPS"));
     88    ppStackMemDump("reject");
    21289
    213     float threshold = psMetadataLookupF32(NULL, recipe, "THRESHOLD.MASK"); // Threshold for mask deconvolution
    214     float imageRej = psMetadataLookupF32(NULL, recipe, "IMAGE.REJ"); // Maximum fraction of image to reject
    215                                                                      // before rejecting entire image
    216     float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"
    217 
    218     const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); // Filename for statistics
    219     psMetadata *stats = NULL;           // Container for statistics
    220     FILE *statsFile = NULL;             // File stream for statistics
    221     if (statsName && strlen(statsName) > 0) {
    222         psString resolved = pmConfigConvertFilename(statsName, config, true, true); // Resolved filename
    223         statsFile = fopen(resolved, "w");
    224         if (!statsFile) {
    225             psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
    226             psFree(resolved);
    227             return false;
    228         } else {
    229             stats = psMetadataAlloc();
    230         }
    231         psFree(resolved);
    232     }
    233 
    234     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file
    235     if (!output) {
    236         psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!");
    237         return false;
    238     }
    239     int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
    240 
    241     psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    242     int overlap = 2 * psMetadataLookupS32(NULL, ppsub,
    243                                           "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans
    244 
    245     if (!pmConfigMaskSetBits(NULL, NULL, config)) {
    246         psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value.");
    247         return false;
    248     }
    249 
    250     memDump("start");
    251 
    252     psLogMsg("ppStack", PS_LOG_INFO, "Stage 0 : Initialization and Configuration : %f sec", psTimerClear("PPSTACK_STEPS"));
    253 
    254     // Preparation iteration: Load the sources, and get a target PSF model
    255     psTrace("ppStack", 1, "Determining target PSF....\n");
    256     psArray *sourceLists = psArrayAlloc(num); // Individual lists of sources for matching
    257     pmPSF *targetPSF = NULL;            // Target PSF
    258     float sumExposure = NAN;            // Sum of exposure times
    259     if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
    260         pmFPAfileActivate(config->files, false, NULL);
    261         fileActivation(config, prepareFiles, true);
    262         pmFPAview *view = filesIterateDown(config);
    263         if (!view) {
    264             return false;
    265         }
    266 
    267         // Generate list of PSFs
    268         psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
    269                                                                "^PPSTACK.INPUT$"); // Iterator
    270         psMetadataItem *fileItem; // Item from iteration
    271         psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
    272         int numCols = 0, numRows = 0;   // Size of image
    273         int index = 0;                  // Index for file
    274         while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
    275             assert(fileItem->type == PS_DATA_UNKNOWN);
    276             pmFPAfile *inputFile = fileItem->data.V; // An input file
    277             pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
    278             pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
    279             if (!psf) {
    280                 psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
    281                 psFree(view);
    282                 psFree(sourceLists);
    283                 psFree(fileIter);
    284                 psFree(psfs);
    285                 return false;
    286             }
    287             psfs->data[index] = psMemIncrRefCounter(psf);
    288             psMetadataRemoveKey(chip->analysis, "PSPHOT.PSF");
    289 
    290             pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
    291             pmHDU *hdu = pmHDUFromCell(cell);
    292             assert(hdu && hdu->header);
    293             int naxis1 = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns
    294             int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows
    295             if (naxis1 <= 0 || naxis2 <= 0) {
    296                 psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF.");
    297                 psFree(view);
    298                 psFree(sourceLists);
    299                 psFree(fileIter);
    300                 psFree(psfs);
    301                 return false;
    302             }
    303             if (numCols == 0 && numRows == 0) {
    304                 numCols = naxis1;
    305                 numRows = naxis2;
    306             }
    307 
    308             pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
    309             psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
    310             if (!sources) {
    311                 psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
    312                 return NULL;
    313             }
    314             sourceLists->data[index] = psMemIncrRefCounter(sources);
    315 
    316             // Re-do photometry if we don't trust the source lists
    317             if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
    318                 psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
    319                 pmFPAfileActivate(config->files, false, NULL);
    320                 fileActivationSingle(config, convolveFiles, true, index);
    321                 pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
    322                 pmFPAview *view = filesIterateDown(config);
    323                 if (!view) {
    324                     psFree(sourceLists);
    325                     psFree(targetPSF);
    326                     return false;
    327                 }
    328 
    329                 pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
    330 
    331                 if (!ppStackInputPhotometry(ro, sources, config)) {
    332                     psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
    333                     psFree(sourceLists);
    334                     psFree(targetPSF);
    335                     return false;
    336                 }
    337 
    338                 psFree(view);
    339                 if (!filesIterateUp(config)) {
    340                     psFree(sourceLists);
    341                     psFree(targetPSF);
    342                     return false;
    343                 }
    344                 pmFPAfileActivate(config->files, false, NULL);
    345                 fileActivation(config, prepareFiles, true);
    346             }
    347 
    348             index++;
    349         }
    350         psFree(fileIter);
    351 
    352         // Zero point calibration
    353         sumExposure = ppStackSourcesTransparency(sourceLists, view, config);
    354         if (!isfinite(sumExposure) || sumExposure <= 0) {
    355             psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
    356             psFree(sourceLists);
    357             psFree(targetPSF);
    358             return false;
    359         }
    360 
    361         // Generate target PSF
    362         targetPSF = ppStackPSF(config, numCols, numRows, psfs);
    363         psFree(psfs);
    364         if (!targetPSF) {
    365             psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
    366             psFree(sourceLists);
    367             psFree(view);
    368             return false;
    369         }
    370         psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN,
    371                          "Target PSF for stack", targetPSF);
    372 
    373         pmChip *outChip = pmFPAviewThisChip(view, output->fpa); // Output chip
    374         psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,
    375                          "Target PSF", targetPSF);
    376         outChip->data_exists = true;
    377 
    378         psFree(view);
    379         if (!filesIterateUp(config)) {
    380             return false;
    381         }
    382     }
    383 
    384     psLogMsg("ppStack", PS_LOG_INFO, "Stage 1 : Load Sources and Generate Target PSF : %f sec", psTimerClear("PPSTACK_STEPS"));
    385     memDump("psf");
    386 
    387     psArray *imageNames = psArrayAlloc(num);
    388     psArray *maskNames = psArrayAlloc(num);
    389     psArray *varianceNames = psArrayAlloc(num);
    390     for (int i = 0; i < num; i++) {
    391         psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
    392         psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage);
    393         psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask);
    394         psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance);
    395         psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName);
    396         imageNames->data[i] = imageName;
    397         maskNames->data[i] = maskName;
    398         varianceNames->data[i] = varianceName;
    399     }
    400     // Free the outputName that we grabbed above to set up tempName.
    401     psFree(outputName);
    402 
    403     // Generate convolutions and write them to disk
    404     psTrace("ppStack", 1, "Convolving inputs to target PSF....\n");
    405     psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again
    406     psArray *subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking
    407     psArray *subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking
    408     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
    409     int numGood = 0;                    // Number of good frames
    410     int numCols = 0, numRows = 0;       // Size of image
    411     psVector *inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs
    412     psVectorInit(inputMask, 0);
    413     psVector *matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching
    414     psVectorInit(matchChi2, NAN);
    415     psVector *weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
    416     psVectorInit(weightings, NAN);
    417     psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
    418     psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
    419     psArray *covariances = psArrayAlloc(num); // Covariance matrices
    420     for (int i = 0; i < num; i++) {
    421         psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num);
    422         pmFPAfileActivate(config->files, false, NULL);
    423         fileActivationSingle(config, convolveFiles, true, i);
    424         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
    425         pmFPAview *view = filesIterateDown(config);
    426         if (!view) {
    427             psFree(sourceLists);
    428             psFree(targetPSF);
    429             psFree(rng);
    430             psFree(inputMask);
    431             psFree(matchChi2);
    432             psFree(fpaList);
    433             psFree(cellList);
    434             psFree(covariances);
    435             return false;
    436         }
    437 
    438         pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Input readout
    439         psFree(view);
    440 
    441         if (numCols == 0 && numRows == 0) {
    442             numCols = readout->image->numCols;
    443             numRows = readout->image->numRows;
    444         } else if (numCols != readout->image->numCols || numRows != readout->image->numRows) {
    445             psError(PS_ERR_UNKNOWN, true, "Sizes of input images don't match: %dx%d vs %dx%d",
    446                     readout->image->numCols, readout->image->numRows, numCols, numRows);
    447             psFree(sourceLists);
    448             psFree(targetPSF);
    449             psFree(rng);
    450             psFree(inputMask);
    451             psFree(matchChi2);
    452             psFree(fpaList);
    453             psFree(cellList);
    454             psFree(covariances);
    455             return false;
    456         }
    457 
    458         // Background subtraction, scaling and normalisation is performed automatically by the image matching
    459         psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
    460         psTimerStart("PPSTACK_MATCH");
    461 
    462         if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i], &weightings->data.F32[i],
    463                           sourceLists->data[i], targetPSF, rng, config)) {
    464             psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
    465             inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_MATCH;
    466             psErrorClear();
    467             continue;
    468         }
    469         covariances->data[i] = psMemIncrRefCounter(readout->covariance);
    470 
    471         if (stats) {
    472             pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
    473                                                                 PM_SUBTRACTION_ANALYSIS_KERNEL);// Conv kernel
    474             psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", PS_META_DUPLICATE_OK,
    475                              "Time to match PSF", psTimerMark("PPSTACK_MATCH"));
    476             psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK,
    477                              "Mean deviation for stamps", kernels->mean);
    478             psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK,
    479                              "RMS deviation for stamps", kernels->rms);
    480             psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK,
    481                              "Number of stamps", kernels->numStamps);
    482             float deconv = psMetadataLookupF32(NULL, readout->analysis,
    483                                                PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Deconvolution fraction
    484             psMetadataAddF32(stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK,
    485                              "Deconvolution fraction for kernel", deconv);
    486             psMetadataAddF32(stats, PS_LIST_TAIL, "PPSTACK.WEIGHTING", PS_META_DUPLICATE_OK,
    487                              "Weighting for image", weightings->data.F32[i]);
    488         }
    489         psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH"));
    490 
    491         subRegions->data[i] = regions;
    492         subKernels->data[i] = kernels;
    493 
    494         // Write the temporary convolved files
    495         pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
    496         assert(hdu);
    497         writeImage(imageNames->data[i], hdu->header, readout->image, config);
    498         psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
    499         pmConfigMaskWriteHeader(config, maskHeader);
    500         writeImage(maskNames->data[i], maskHeader, readout->mask, config);
    501         psFree(maskHeader);
    502         psImageCovarianceTransfer(readout->variance, readout->covariance);
    503         writeImage(varianceNames->data[i], hdu->header, readout->variance, config);
    504 #ifdef TESTING
    505         {
    506             psString name = NULL;
    507             psStringAppend(&name, "covariance_%d.fits", i);
    508             writeImage(name, hdu->header, readout->covariance->image, config);
    509             psFree(name);
    510         }
    511 #endif
    512 
    513         pmCell *inCell = readout->parent; // Input cell
    514 
    515         psListAdd(cellList, PS_LIST_TAIL, inCell);
    516         psListAdd(fpaList, PS_LIST_TAIL, inCell->parent->parent);
    517 
    518         cells->data[i] = psMemIncrRefCounter(inCell);
    519         if (!filesIterateUp(config)) {
    520             psFree(fpaList);
    521             psFree(cellList);
    522             psFree(covariances);
    523             return false;
    524         }
    525         numGood++;
    526 
    527         memDump("match");
    528     }
    529     psFree(sourceLists);
    530     psFree(targetPSF);
    531     psFree(rng);
    532 
    533     psLogMsg("ppStack", PS_LOG_INFO, "Stage 2 : Generate Convolutions and Save : %f sec", psTimerClear("PPSTACK_STEPS"));
    534 
    535     if (numGood == 0) {
    536         psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
    537         psFree(subKernels);
    538         psFree(subRegions);
    539         psFree(cells);
    540         psFree(inputMask);
    541         psFree(matchChi2);
    542         psFree(fpaList);
    543         psFree(cellList);
    544         psFree(covariances);
    545         return false;
    546     }
    547 
    548 
    549     // Update concepts for output
    550     {
    551         pmFPAview view;                 // View for output
    552         view.chip = view.cell = view.readout = 0;
    553         pmCell *outCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.OUTPUT"); // Output cell
    554         pmFPA *outFPA = outCell->parent->parent; // Output FPA
    555         pmConceptsAverageFPAs(outFPA, fpaList);
    556         pmConceptsAverageCells(outCell, cellList, NULL, NULL, true);
    557         psFree(fpaList);
    558         psFree(cellList);
    559 
    560         // Update the value of a concept
    561 #define UPDATE_CONCEPT(SOURCE, NAME, VALUE) { \
    562             psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \
    563             psAssert(item, "Concept should be present"); \
    564             psAssert(item->type == PS_TYPE_F32, "Concept should be F32"); \
    565             item->data.F32 = VALUE; \
    566         }
    567 
    568         UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",  sumExposure);
    569         UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", sumExposure);
    570         UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
    571     }
    572 
    573 
    574     // Reject images out-of-hand on the basis of their match chi^2
    575     {
    576         psVector *values = psVectorAllocEmpty(num, PS_TYPE_F32); // Values to sort
    577         for (int i = 0; i < num; i++) {
    578             if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
    579                 continue;
    580             }
    581             values->data.F32[values->n++] = matchChi2->data.F32[i];
    582         }
    583         assert(values->n == numGood);
    584         if (!psVectorSortInPlace(values)) {
    585             psError(PS_ERR_UNKNOWN, false, "Unable to sort vector.");
    586             psFree(subKernels);
    587             psFree(subRegions);
    588             psFree(cells);
    589             psFree(inputMask);
    590             psFree(matchChi2);
    591             psFree(values);
    592             psFree(covariances);
    593             return false;
    594         }
    595         float median = numGood % 2 ? values->data.F32[numGood / 2] :
    596             0.5 * (values->data.F32[numGood / 2 - 1] + values->data.F32[numGood / 2]);
    597 
    598         float rms = 0.74 * (values->data.F32[numGood * 3 / 4] -
    599                             values->data.F32[numGood / 4]); // Estimated RMS from interquartile range
    600 
    601         psFree(values);
    602 
    603         float rej = psMetadataLookupF32(NULL, recipe, "MATCH.REJ"); // Rejection threshold (stdevs)
    604         if (isfinite(rej)) {
    605             float thresh = median + rej * rms; // Threshold for rejection
    606             psLogMsg("ppStack", PS_LOG_INFO, "chi^2 rejection threshold = %f + %f * %f = %f",
    607                      median, rej, rms, thresh);
    608 
    609             int numRej = 0;                 // Number rejected
    610             numGood = 0;                    // Number of good images
    611             for (int i = 0; i < num; i++) {
    612               if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
    613                     continue;
    614                 }
    615                 if (matchChi2->data.F32[i] > thresh) {
    616                     numRej++;
    617                     inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CHI2;
    618                     psLogMsg("ppStack", PS_LOG_INFO, "Rejecting image %d because of large matching chi^2: %f",
    619                              i, matchChi2->data.F32[i]);
    620                 } else {
    621                     psLogMsg("ppStack", PS_LOG_INFO, "Image %d has matching chi^2: %f",
    622                              i, matchChi2->data.F32[i]);
    623                     numGood++;
    624                 }
    625             }
    626         }
    627     }
    628 
    629     if (numGood == 0) {
    630         psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
    631         psFree(subKernels);
    632         psFree(subRegions);
    633         psFree(cells);
    634         psFree(inputMask);
    635         psFree(matchChi2);
    636         psFree(covariances);
    637         return false;
    638     }
    639 
    640 #ifdef TESTING
    641     psTraceSetLevel("psModules.imcombine", 7);
    642 #endif
    643 
    644     psLogMsg("ppStack", PS_LOG_INFO, "Stage 3 : Basic Image Rejection  : %f sec", psTimerClear("PPSTACK_STEPS"));
    645 
    646     // Start threading
    647     ppStackThreadInit();
    648     ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, varianceNames,
    649                                                       covariances, config);
    650     if (!stack) {
    651         psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
    652         psFree(subKernels);
    653         psFree(subRegions);
    654         psFree(inputMask);
    655         psFree(matchChi2);
    656         psFree(cells);
    657         psFree(covariances);
    658         return false;
    659     }
    660 
    661     psTimerStart("PPSTACK_INITIAL");
    662 
    663     memDump("preinitial");
    664 
    665     // Stack the convolved files
    666     psTrace("ppStack", 1, "Initial stack of convolved images....\n");
    667     pmReadout *outRO = NULL;            // Output readout
    668     pmFPAview *view = NULL;             // View to readout
    669     psArray *inspect = NULL;            // Array of arrays of pixels to inspect
    670     {
    671         int row0, col0;                 // Offset for readout
    672         int numCols, numRows;           // Size of readout
    673         ppStackThread *thread = stack->threads->data[0]; // Representative thread
    674         if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) {
    675             psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
    676             psFree(subKernels);
    677             psFree(subRegions);
    678             psFree(stack);
    679             psFree(inputMask);
    680             psFree(matchChi2);
    681             psFree(cells);
    682             psFree(covariances);
    683             return false;
    684         }
    685 
    686         pmFPAfileActivate(config->files, false, NULL);
    687         fileActivation(config, combineFiles, true);
    688         view = filesIterateDown(config);
    689         if (!view) {
    690             psFree(subKernels);
    691             psFree(subRegions);
    692             psFree(stack);
    693             psFree(inputMask);
    694             psFree(matchChi2);
    695             psFree(cells);
    696             psFree(covariances);
    697             return false;
    698         }
    699 
    700         pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
    701         outRO = pmReadoutAlloc(outCell); // Output readout
    702 
    703         psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    704         psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    705         if (!pmReadoutStackDefineOutput(outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
    706             psError(PS_ERR_UNKNOWN, false, "Unable to prepare output.");
    707             psFree(subKernels);
    708             psFree(subRegions);
    709             psFree(stack);
    710             psFree(inputMask);
    711             psFree(matchChi2);
    712             psFree(view);
    713             psFree(outRO);
    714             psFree(cells);
    715             psFree(covariances);
    716             return false;
    717         }
    718 
    719         psFree(cells);
    720 
    721         bool status;                    // Status of read
    722         int numChunk;                   // Number of chunks
    723         for (numChunk = 0; true; numChunk++) {
    724             ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, overlap);
    725             if (!status) {
    726                 // Something went wrong
    727                 psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
    728                 psFree(subKernels);
    729                 psFree(subRegions);
    730                 psFree(stack);
    731                 psFree(inputMask);
    732                 psFree(matchChi2);
    733                 psFree(view);
    734                 psFree(outRO);
    735                 psFree(covariances);
    736                 return false;
    737             }
    738             if (!thread) {
    739                 // Nothing more to read
    740                 break;
    741             }
    742 
    743             // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)
    744             psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start
    745             psArrayAdd(job->args, 1, thread);
    746             psArrayAdd(job->args, 1, config);
    747             psArrayAdd(job->args, 1, outRO);
    748             psArrayAdd(job->args, 1, inputMask);
    749             psArrayAdd(job->args, 1, weightings);
    750             psArrayAdd(job->args, 1, matchChi2);
    751             if (!psThreadJobAddPending(job)) {
    752                 psFree(job);
    753                 psFree(subKernels);
    754                 psFree(subRegions);
    755                 psFree(stack);
    756                 psFree(inputMask);
    757                 psFree(weightings);
    758                 psFree(matchChi2);
    759                 psFree(view);
    760                 psFree(outRO);
    761                 psFree(covariances);
    762                 return false;
    763             }
    764             psFree(job);
    765         }
    766 
    767         if (!psThreadPoolWait(false)) {
    768             psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination.");
    769             psFree(subKernels);
    770             psFree(subRegions);
    771             psFree(stack);
    772             psFree(inputMask);
    773             psFree(matchChi2);
    774             psFree(view);
    775             psFree(outRO);
    776             psFree(covariances);
    777             return false;
    778         }
    779 
    780         // Harvest the jobs, gathering the inspection lists
    781         inspect = psArrayAlloc(num);
    782         for (int i = 0; i < num; i++) {
    783             if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    784                 continue;
    785             }
    786             inspect->data[i] = psArrayAllocEmpty(numChunk);
    787         }
    788         psThreadJob *job;               // Completed job
    789         while ((job = psThreadJobGetDone())) {
    790             psAssert(strcmp(job->type, "PPSTACK_INITIAL_COMBINE") == 0,
    791                      "Job has incorrect type: %s", job->type);
    792             psArray *results = job->results; // Results of job
    793             for (int i = 0; i < num; i++) {
    794                 if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    795                     continue;
    796                 }
    797                 inspect->data[i] = psArrayAdd(inspect->data[i], 1, results->data[i]);
    798             }
    799             psFree(job);
    800         }
    801 
    802         memDump("initial");
    803     }
    804 
    805 #ifdef TESTING
    806     writeImage("combined_initial.fits", NULL, outRO->image, config);
    807 #endif
    808 
    809     if (stats) {
    810         psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_INITIAL", 0,
    811                          "Time to make initial stack", psTimerMark("PPSTACK_INITIAL"));
    812     }
    813     psLogMsg("ppStack", PS_LOG_INFO, "Stage 4 : Make Initial Stack : %f sec", psTimerClear("PPSTACK_STEPS"));
    814     psLogMsg("ppStack", PS_LOG_INFO, "Time to make initial stack: %f sec", psTimerClear("PPSTACK_INITIAL"));
    815 
    816     psTimerStart("PPSTACK_REJECT");
    817 
    818     // Pixel rejection
    819     psArray *rejected = psArrayAlloc(num);
    820     {
    821         int numRejected = 0;        // Number of inputs rejected completely
    822 
    823         // Count images rejected out of hand
    824         for (int i = 0; i < num; i++) {
    825             if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    826                 numRejected++;
    827             }
    828         }
    829 
    830         for (int i = 0; i < num; i++) {
    831             if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    832                 continue;
    833             }
    834 
    835             psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start
    836             psArrayAdd(job->args, 1, inspect);
    837             PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
    838             if (!psThreadJobAddPending(job)) {
    839                 psFree(job);
    840                 psFree(subKernels);
    841                 psFree(subRegions);
    842                 psFree(stack);
    843                 psFree(inputMask);
    844                 psFree(view);
    845                 psFree(outRO);
    846                 psFree(inspect);
    847                 psFree(rejected);
    848                 psFree(covariances);
    849                 psFree(matchChi2);
    850                 return false;
    851             }
    852             psFree(job);
    853         }
    854 
    855         if (!psThreadPoolWait(true)) {
    856             psError(PS_ERR_UNKNOWN, false, "Unable to concatenate inspection lists.");
    857             psFree(subKernels);
    858             psFree(subRegions);
    859             psFree(stack);
    860             psFree(inputMask);
    861             psFree(view);
    862             psFree(outRO);
    863             psFree(inspect);
    864             psFree(rejected);
    865             psFree(covariances);
    866             psFree(matchChi2);
    867             return false;
    868         }
    869 
    870         if (psMetadataLookupS32(NULL, config->arguments, "-threads") > 0) {
    871             pmStackRejectThreadsInit();
    872         }
    873 
    874         int stride = psMetadataLookupS32(NULL, ppsub, "STRIDE"); // Size of convolution patches
    875 
    876         // Reject bad pixels
    877         for (int i = 0; i < num; i++) {
    878             if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    879                 continue;
    880             }
    881             psTimerStart("PPSTACK_REJECT");
    882 
    883 #ifdef TESTING
    884             {
    885                 psImage *mask = psPixelsToMask(NULL, inspect->data[i],
    886                                                psRegionSet(0, numCols - 1, 0, numRows - 1),
    887                                                0xff); // Mask image
    888                 psString name = NULL;           // Name of image
    889                 psStringAppend(&name, "inspect_%03d.fits", i);
    890                 psFits *fits = psFitsOpen(name, "w");
    891                 psFree(name);
    892                 psFitsWriteImage(fits, NULL, mask, 0, NULL);
    893                 psFree(mask);
    894                 psFitsClose(fits);
    895             }
    896 #endif
    897 
    898             psPixels *reject = pmStackReject(inspect->data[i], numCols, numRows, threshold, poorFrac, stride,
    899                                              subRegions->data[i], subKernels->data[i]); // Rejected pixels
    900 
    901 #ifdef TESTING
    902             {
    903                 psImage *mask = psPixelsToMask(NULL, reject, psRegionSet(0, numCols - 1, 0, numRows - 1),
    904                                                0xff); // Mask image
    905                 psString name = NULL;           // Name of image
    906                 psStringAppend(&name, "reject_%03d.fits", i);
    907                 psFits *fits = psFitsOpen(name, "w");
    908                 psFree(name);
    909                 psFitsWriteImage(fits, NULL, mask, 0, NULL);
    910                 psFree(mask);
    911                 psFitsClose(fits);
    912             }
    913 #endif
    914 
    915             psFree(inspect->data[i]);
    916             inspect->data[i] = NULL;
    917 
    918             if (!reject) {
    919                 psWarning("Rejection on image %d didn't work --- reject entire image.", i);
    920                 numRejected++;
    921                 inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_REJECT;
    922             } else {
    923                 float frac = reject->n / (float)(numCols * numRows); // Pixel fraction
    924                 psLogMsg("ppStack", PS_LOG_INFO, "%ld pixels rejected from image %d (%.1f%%)",
    925                          reject->n, i, frac * 100.0);
    926                 if (frac > imageRej) {
    927                     psWarning("Image %d rejected completely because rejection fraction (%.3f) "
    928                               "exceeds limit (%.3f)", i, frac, imageRej);
    929                     psFree(reject);
    930                     // reject == NULL means reject image completely
    931                     reject = NULL;
    932                     inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_BAD;
    933                     numRejected++;
    934                 }
    935             }
    936 
    937             // Images without a list of rejected pixels (the list may be empty) are rejected completely
    938             rejected->data[i] = reject;
    939 
    940             if (stats) {
    941                 psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_REJECT", PS_META_DUPLICATE_OK,
    942                                  "Time to perform rejection", psTimerMark("PPSTACK_REJECT"));
    943                 psMetadataAddS32(stats, PS_LIST_TAIL, "REJECT_PIXELS", PS_META_DUPLICATE_OK,
    944                                  "Number of pixels rejected", reject ? reject->n : 0);
    945             }
    946             psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i,
    947                      psTimerClear("PPSTACK_REJECT"));
    948         }
    949 
    950         psFree(inspect);
    951         psFree(subKernels);
    952         psFree(subRegions);
    953 
    954         if (numRejected >= num - 1) {
    955             psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed.");
    956             psFree(stack);
    957             psFree(rejected);
    958             psFree(inputMask);
    959             psFree(view);
    960             psFree(outRO);
    961             psFree(covariances);
    962             psFree(matchChi2);
    963             return false;
    964         }
    965 
    966         if (stats) {
    967             psMetadataAddS32(stats, PS_LIST_TAIL, "REJECT_IMAGES", 0,
    968                              "Number of images rejected completely", numRejected);
    969         }
    970     }
    971 
    972     psLogMsg("ppStack", PS_LOG_INFO, "Stage 5 : Pixel Rejection : %f sec", psTimerClear("PPSTACK_STEPS"));
    973     psTimerStart("PPSTACK_FINAL");
    974 
    975     memDump("reject");
    97690
    97791    // Final combination
    97892    psTrace("ppStack", 2, "Final stack of convolved images....\n");
    979     {
    980         stack->lastScan = 0;            // Reset read
    981         bool status;                    // Status of read
    982         for (int numChunk = 0; true; numChunk++) {
    983             ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, 0);
    984             if (!status) {
    985                 // Something went wrong
    986                 psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
    987                 psFree(stack);
    988                 psFree(rejected);
    989                 psFree(inputMask);
    990                 psFree(view);
    991                 psFree(outRO);
    992                 psFree(covariances);
    993                 psFree(matchChi2);
    994                 return false;
    995             }
    996             if (!thread) {
    997                 // Nothing more to read
    998                 break;
    999             }
     93    if (!ppStackCombineFinal(stack, options, config)) {
     94        psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination.");
     95        psFree(stack);
     96        psFree(options);
     97        return false;
     98    }
     99    psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
     100    ppStackMemDump("final");
    1000101
    1001             // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
    1002             psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
    1003             psArrayAdd(job->args, 1, thread);
    1004             psArrayAdd(job->args, 1, config);
    1005             psArrayAdd(job->args, 1, outRO);
    1006             psArrayAdd(job->args, 1, inputMask);
    1007             psArrayAdd(job->args, 1, rejected);
    1008             psArrayAdd(job->args, 1, weightings);
    1009             psArrayAdd(job->args, 1, matchChi2);
    1010             if (!psThreadJobAddPending(job)) {
    1011                 psFree(job);
    1012                 psFree(stack);
    1013                 psFree(rejected);
    1014                 psFree(inputMask);
    1015                 psFree(view);
    1016                 psFree(outRO);
    1017                 psFree(covariances);
    1018                 psFree(matchChi2);
    1019                 return false;
    1020             }
    1021             psFree(job);
    1022         }
    1023102
    1024         if (!psThreadPoolWait(true)) {
    1025             psError(PS_ERR_UNKNOWN, false, "Unable to do final combination.");
    1026             psFree(stack);
    1027             psFree(inputMask);
    1028             psFree(rejected);
    1029             psFree(view);
    1030             psFree(outRO);
    1031             psFree(covariances);
    1032             psFree(matchChi2);
    1033             return false;
    1034         }
     103    // Clean up
     104    psTrace("ppStack", 2, "Cleaning up after combination....\n");
     105    if (!ppStackCombineFinal(stack, options, config)) {
     106        psError(PS_ERR_UNKNOWN, false, "Unable to clean up.");
     107        psFree(stack);
     108        psFree(options);
     109        return false;
    1035110    }
     111    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
     112    ppStackMemDump("cleanup");
    1036113
    1037     memDump("final");
    1038 
    1039 #ifdef TESTING
    1040     writeImage("combined_final.fits", NULL, outRO->image, config);
    1041 #endif
    1042 
    1043     // Sum covariance matrices
    1044     for (int i = 0; i < num; i++) {
    1045         if (inputMask->data.U8[i]) {
    1046             psFree(covariances->data[i]);
    1047             covariances->data[i] = NULL;
    1048         }
    1049     }
    1050     outRO->covariance = psImageCovarianceSum(covariances);
    1051     psFree(covariances);
    1052     psFree(matchChi2);
    1053 
    1054     if (stats) {
    1055         psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_FINAL", 0,
    1056                          "Time to make final stack", psTimerMark("PPSTACK_FINAL"));
    1057     }
    1058 
    1059     psLogMsg("ppStack", PS_LOG_INFO, "Stage 6 : Final Stack : %f sec", psTimerClear("PPSTACK_STEPS"));
    1060     psLogMsg("ppStack", PS_LOG_INFO, "Time to make final stack: %f sec", psTimerClear("PPSTACK_FINAL"));
    1061 
    1062     psTrace("ppStack", 2, "Cleaning up after combination....\n");
    1063 
    1064 #if 0
    1065     // Ensure masked regions really look masked
    1066     {
    1067         psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
    1068         psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    1069         if (!pmReadoutMaskApply(outRO, maskBad)) {
    1070             psWarning("Unable to apply mask");
    1071         }
    1072     }
    1073 #endif
    1074 
    1075     // Close up
    1076     bool wcsDone = false;           // Have we done the WCS?
    1077     for (int i = 0; i < num; i++) {
    1078         if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    1079             continue;
    1080         }
    1081 
    1082         ppStackThread *thread = stack->threads->data[0]; // Representative stack
    1083         pmReadout *inRO = thread->readouts->data[i]; // Template readout
    1084         if (inRO && !wcsDone) {
    1085             // Copy astrometry over
    1086             wcsDone = true;
    1087             pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
    1088             pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU
    1089             pmChip *outChip = outRO->parent->parent; // Output chip
    1090             pmFPA *outFPA = outChip->parent; // Output FPA
    1091             if (!outHDU || !inHDU) {
    1092                 psWarning("Unable to find HDU at FPA level to copy astrometry.");
    1093             } else {
    1094                 if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
    1095                     psErrorClear();
    1096                     psWarning("Unable to read WCS astrometry from input FPA.");
    1097                     wcsDone = false;
    1098                 } else {
    1099                     if (!outHDU->header) {
    1100                         outHDU->header = psMetadataAlloc();
    1101                     }
    1102                     if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
    1103                         psErrorClear();
    1104                         psWarning("Unable to write WCS astrometry to output FPA.");
    1105                         wcsDone = false;
    1106                     }
    1107                 }
    1108             }
    1109         }
    1110 
    1111         if (tempDelete) {
    1112             psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
    1113             psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
    1114             psString varianceResolved = pmConfigConvertFilename(varianceNames->data[i], config, false, false);
    1115             if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||
    1116                 unlink(varianceResolved) == -1) {
    1117                 psWarning("Unable to delete temporary files for image %d", i);
    1118             }
    1119             psFree(imageResolved);
    1120             psFree(maskResolved);
    1121             psFree(varianceResolved);
    1122         }
    1123     }
    1124     psFree(imageNames);
    1125     psFree(maskNames);
    1126     psFree(varianceNames);
    1127 
    1128     psFree(inputMask);
    1129114    psFree(stack);
    1130115
    1131     memDump("cleanup");
    1132116
    1133     // Generate binned JPEGs
    1134     {
    1135         int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level
    1136         int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level
    1137 
    1138         // Target cells
    1139         pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG1");
    1140         pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG2");
    1141         psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    1142 
    1143         pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts
    1144         if (!pmReadoutRebin(ro1, outRO, maskValue, bin1, bin1) || !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) {
    1145             psError(PS_ERR_UNKNOWN, false, "Unable to bin output.");
    1146             psFree(ro1);
    1147             psFree(ro2);
    1148             psFree(outRO);
    1149             return false;
    1150         }
    1151         psFree(ro1);
    1152         psFree(ro2);
    1153     }
    1154 
    1155     psLogMsg("ppStack", PS_LOG_INFO, "Stage 7 : WCS & JPEGS : %f sec", psTimerClear("PPSTACK_STEPS"));
    1156 
    1157     if (psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
    1158         psTrace("ppStack", 1, "Photometering stacked image....\n");
    1159 
    1160         psTimerStart("PPSTACK_PHOT");
    1161 
    1162         fileActivation(config, combineFiles, false);
    1163         fileActivation(config, photFiles, true);
    1164         pmFPAview *photView = filesIterateDown(config);
    1165         if (!ppStackPhotometry(config, outRO, photView)) {
    1166             psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
    1167             psFree(outRO);
    1168             psFree(photView);
    1169             return false;
    1170         }
    1171         psFree(photView);
    1172 
    1173         fileActivation(config, combineFiles, true);
    1174 
    1175         if (stats) {
    1176             pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File
    1177             pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
    1178             psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
    1179             psMetadataAddS32(stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", sources->n);
    1180             psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_PHOT", 0,
    1181                              "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
    1182         }
    1183         psLogMsg("ppStack", PS_LOG_INFO, "Time to do photometry: %f sec", psTimerClear("PPSTACK_PHOT"));
    1184     }
    1185 
    1186     psLogMsg("ppStack", PS_LOG_INFO, "Stage 8 : Photometry Analysis : %f sec", psTimerClear("PPSTACK_STEPS"));
    1187 
    1188     psThreadPoolFinalize();
    1189 
    1190     memDump("phot");
    1191 
    1192     // Statistics on output
    1193     if (stats) {
    1194         psTrace("ppStack", 1, "Gathering statistics on stacked image....\n");
    1195         psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
    1196         psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    1197 
    1198         ppStatsFPA(stats, outRO->parent->parent->parent, view, maskBad, config);
    1199     }
    1200 
    1201     memDump("stats");
    1202 
    1203     // Put version information into the header
    1204     pmHDU *hdu = pmHDUFromCell(outRO->parent);
    1205     if (!hdu) {
    1206         psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
    1207         psFree(outRO);
    1208         psFree(view);
     117    // Photometry
     118    psTrace("ppStack", 1, "Photometering stacked image....\n");
     119    if (!ppStackPhotometry(options, config)) {
     120        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry.");
     121        psFree(options);
    1209122        return false;
    1210123    }
    1211     if (!hdu->header) {
    1212         hdu->header = psMetadataAlloc();
    1213     }
    1214     ppStackVersionMetadata(hdu->header);
     124    psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
     125    ppStackMemDump("photometry");
    1215126
    1216     psFree(outRO);
    1217     psFree(view);
    1218127
    1219     // Write out summary statistics
    1220     if (stats) {
    1221         psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_STACK", 0,
    1222                          "Time to do photometry", psTimerClear("PPSTACK_TOTAL"));
    1223 
    1224         const char *statsMDC = psMetadataConfigFormat(stats);
    1225         if (!statsMDC || strlen(statsMDC) == 0) {
    1226             psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
    1227         } else {
    1228             fprintf(statsFile, "%s", statsMDC);
    1229         }
    1230         psFree((void *)statsMDC);
    1231         fclose(statsFile);
    1232 
    1233         psFree(stats);
    1234     }
    1235 
    1236     // Write out the output files
    1237     if (!filesIterateUp(config)) {
     128    // Finish up
     129    psTrace("ppStack", 1, "Finishing up....\n");
     130    if (!ppStackFinish(options, config)) {
     131        psError(PS_ERR_UNKNOWN, false, "Unable to finish up.");
     132        psFree(options);
    1238133        return false;
    1239134    }
     135    psLogMsg("ppStack", PS_LOG_INFO, "Stage 9: Final output: %f sec", psTimerClear("PPSTACK_STEPS"));
     136    ppStackMemDump("finish");
    1240137
    1241     psLogMsg("ppStack", PS_LOG_INFO, "Stage 9 : Final Output : %f sec", psTimerClear("PPSTACK_STEPS"));
    1242 
    1243     memDump("finish");
    1244 
     138    psFree(options);
    1245139    return true;
    1246140}

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

    r21525 r23352  
    145145
    146146    psImage *binned = psphotBackgroundModel(ro, config); // Binned background model
    147     psImageBinning *binning = psMetadataLookupPtr(NULL, psphotRecipe,
     147    psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis,
    148148                                                  "PSPHOT.BACKGROUND.BINNING"); // Binning for model
    149149    psAssert(binning, "Need binning parameters");
    … …  
    358358                psString name = NULL;
    359359                psStringAppend(&name, "fake_%03d.fits", numInput);
     360                pmStackVisualPlotTestImage(fake->image, name);
    360361                psFits *fits = psFitsOpen(name, "w");
    361362                psFree(name);
    … …  
    367368                psString name = NULL;
    368369                psStringAppend(&name, "real_%03d.fits", numInput);
     370                pmStackVisualPlotTestImage(readout->image, name);
    369371                psFits *fits = psFitsOpen(name, "w");
    370372                psFree(name);
    … …  
    397399                psString name = NULL;
    398400                psStringAppend(&name, "conv_%03d.fits", numInput);
     401                pmStackVisualPlotTestImage(output->image, name);
    399402                psFits *fits = psFitsOpen(name, "w");
    400403                psFree(name);
    … …  
    406409                psString name = NULL;
    407410                psStringAppend(&name, "diff_%03d.fits", numInput);
     411                pmStackVisualPlotTestImage(fake->image, name);
    408412                psFits *fits = psFitsOpen(name, "w");
    409413                psFree(name);
    … …  
    633637        psString name = NULL;
    634638        psStringAppend(&name, "convolved_%03d.fits", numInput);
     639        pmStackVisualPlotTestImage(output->image, name);
    635640        psFits *fits = psFitsOpen(name, "w");
    636641        psFree(name);

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

    r18918 r23352  
    99#include "ppStack.h"
    1010
    11 pmPSF *ppStackPSF(const pmConfig *config, int numCols, int numRows, const psArray *psfs)
     11pmPSF *ppStackPSF(const pmConfig *config, int numCols, int numRows,
     12                  const psArray *psfs, const psVector *inputMask)
    1213{
    1314    // Get the recipe values
    … …  
    1920    const char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
    2021    int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF
     22
     23    for (int i = 0; i < psfs->n; i++) {
     24        if (inputMask->data.U8[i]) {
     25            psFree(psfs->data[i]);
     26            psfs->data[i] = NULL;
     27        }
     28    }
    2129
    2230    // Solve for the target PSF

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

    r21259 r23352  
    99
    1010#include "ppStack.h"
     11#include "ppStackLoop.h"
    1112
    12 #define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
    13                           PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
    14                           PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources
    15 
    16 bool ppStackInputPhotometry(const pmReadout *ro, const psArray *sources, const pmConfig *config)
     13bool ppStackPhotometry(ppStackOptions *options, pmConfig *config)
    1714{
    18     PM_ASSERT_READOUT_NON_NULL(ro, false);
    19     PS_ASSERT_ARRAY_NON_NULL(sources, false);
     15    psAssert(options, "Require options");
     16    psAssert(config, "Require configuration");
    2017
    2118    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    … …  
    2320
    2421    bool mdok;                          // Status of MD lookup
    25 
    26     float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement
    27     if (!mdok) {
    28         psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe");
    29         return false;
    30     }
    31     float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry
    32     if (!mdok) {
    33         psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe");
    34         return false;
    35     }
    36     float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry
    37     if (!mdok) {
    38         psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe");
    39         return false;
     22    if (!psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
     23        // Nothing to do
     24        return true;
    4025    }
    4126
    42     psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in
    43     if (!mdok || !maskValStr) {
    44         psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe");
    45         return false;
    46     }
    47     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask
     27    psTimerStart("PPSTACK_PHOT");
    4828
    49     psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout
     29    ppStackFileActivation(config, PPSTACK_FILES_COMBINE, false);
     30    ppStackFileActivation(config, PPSTACK_FILES_PHOT, true);
     31    pmFPAview *photView = ppStackFilesIterateDown(config); // View to readout
     32    ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true);
    5033
    51     // Measure background
    52     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
    53     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
    54     if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
    55         psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image");
    56         psFree(stats);
    57         psFree(rng);
    58         return false;
    59     }
    60     float bg = stats->robustMedian; // Background level
    61     psFree(stats);
    62     psFree(rng);
     34    pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File for photometry
     35    pmFPACopy(photFile->fpa, options->outRO->parent->parent->parent);
    6336
    64     // Photometer sources
    65     int numCols = image->numCols, numRows = image->numRows; // Size of image
    66     int numSources = sources->n; // Number of sources
    67     int numGood = 0;            // Number of good sources
    68     float maxMag = -INFINITY;   // Maximum magnitude
    69     for (int j = 0; j < numSources; j++) {
    70         pmSource *source = sources->data[j]; // Source of interest
    71         if (source->mode & PHOT_SOURCE_MASK || !isfinite(source->psfMag)) {
    72             source->psfMag = NAN;
    73             continue;
    74         }
    75         float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
    76         int xCentral = x + 0.5, yCentral = y + 0.5; // Central pixel
    77         // Range of integration
    78         int xMin = PS_MAX(0, xCentral - zpRadius), xMax = PS_MIN(numCols - 1, xCentral + zpRadius);
    79         int yMin = PS_MAX(0, yCentral - zpRadius), yMax = PS_MIN(numRows - 1, yCentral + zpRadius);
    80 
    81         // Integrate
    82         double sumImage = 0.0, sumKernel = 0.0; // Sums from integration
    83         int numBadPix = 0;         // Number of bad pixels
    84         for (int v = yMin; v <= yMax; v++) {
    85             float dy2 = PS_SQR(y - v); // Distance from centroid
    86             for (int u = xMin; u <= xMax; u++) {
    87                 if (mask->data.PS_TYPE_IMAGE_MASK_DATA[v][u] & maskVal) {
    88                     numBadPix++;
    89                     continue;
    90                 }
    91                 float dx2 = PS_SQR(x - u); // Distance from centroid
    92                 double kernel = exp(-0.5 * (dx2 + dy2) / PS_SQR(zpSigma)); // Kernel value
    93                 sumImage += (image->data.F32[v][u] - bg) * kernel;
    94                 sumKernel += kernel;
    95             }
    96         }
    97 
    98         if (numBadPix > zpFrac * M_PI * PS_SQR(zpRadius) || !isfinite(sumImage)) {
    99             source->psfMag = NAN;
    100         } else {
    101             source->psfMag = - 2.5 * log10(sumImage / sumKernel * M_PI * PS_SQR(zpRadius));
    102             if (isfinite(source->psfMag)) {
    103                 numGood++;
    104                 maxMag = PS_MAX(maxMag, source->psfMag);
    105             }
    106         }
    107     }
    108     psLogMsg("ppStack", PS_LOG_INFO, "Photometered %d good sources in image; max mag %f", numGood, maxMag);
    109 
    110     return true;
    111 }
    112 
    113 
    114 
    115 bool ppStackPhotometry(pmConfig *config, const pmReadout *readout, const pmFPAview *view)
    116 {
    117     pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
    118     pmFPACopy(photFile->fpa, readout->parent->parent->parent);
    119 
    120     if (psMetadataLookupBool(NULL, config->arguments, "-psphot-visual")) {
    121         psphotSetVisual(true);
     37    if (psMetadataLookupBool(NULL, config->arguments, "-visual")) {
     38        pmVisualSetVisual(true);
    12239    }
    12340
     41    psMetadata *psphot = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
     42
     43#if 0
    12444    // Need to ensure aperture residual is not calculated
    125     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
    12645    psMetadataItem *item = psMetadataLookup(recipe, "MEASURE.APTREND"); // Item determining aptrend
    12746    if (!item) {
    … …  
    13150        item->data.B = false;
    13251    }
     52#endif
    13353
    13454    // set maskValue and markValue in the psphot recipe
    13555    psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    13656    psImageMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking
    137     psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue);
    138     psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, "Bits to use for mark", markValue);
     57    psMetadataAddImageMask(psphot, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE,
     58                            "Bits to mask", maskValue);
     59    psMetadataAddImageMask(psphot, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE,
     60                           "Bits to use for mark", markValue);
    13961
    140     // XXX replace with psphotReadoutKnownSources (config, view)
    141     // XXX this function requires that we construct the input source list and place it on PSPHOT.INPUT.CMF
    142     pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT");
    143     psArray *inSources = psMetadataLookupPtr (NULL, sourcesCell->analysis, "PSPHOT.SOURCES");
     62    pmCell *sourcesCell = pmFPAfileThisCell(config->files, photView, "PPSTACK.OUTPUT");
     63    psArray *inSources = psMetadataLookupPtr(NULL, sourcesCell->analysis, "PSPHOT.SOURCES");
    14464    if (!inSources) {
    14565        psError(PS_ERR_UNKNOWN, false, "Unable to find input sources");
     66        psFree(photView);
    14667        return false;
    14768    }
    14869
    149     if (!psphotReadoutKnownSources(config, view, inSources)) {
     70    if (!psphotReadoutKnownSources(config, photView, inSources)) {
    15071        // Clear the error, so that the output files are written.
    151         psWarning("Unable to perform photometry on stacked image.");
    152         psErrorStackPrint(stderr, "Error stack from photometry:");
    153         psErrorClear();
     72        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on stacked image.");
     73        psFree(photView);
     74        return false;
    15475    }
     76    psFree(photView);
    15577
    15678    if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
    … …  
    16385    pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
    16486
     87    if (options->stats) {
     88        pmReadout *photRO = pmFPAviewThisReadout(photView, photFile->fpa); // Readout with the sources
     89        psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
     90        psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0,
     91                         "Number of sources detected", sources->n);
     92        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0,
     93                         "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
     94    }
     95    psLogMsg("ppStack", PS_LOG_INFO, "Time to do photometry: %f sec", psTimerClear("PPSTACK_PHOT"));
     96
    16597    return true;
    16698}

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

    r21477 r23352  
    88#include <psphot.h>
    99
     10#include "ppStackThread.h"
    1011#include "ppStack.h"
    1112
    … …  
    1617    psArray *args = job->args;          // Arguments
    1718    ppStackThread *thread = args->data[0]; // Thread
    18     pmConfig *config = args->data[1];   // Configuration
    19     pmReadout *outRO = args->data[2];   // Output readout
    20     psVector *mask = args->data[3];     // Mask for inputs
    21     psVector *weightings = args->data[4]; // Weightings (1/noise^2) for each image
    22     psVector *addVariance = args->data[5]; // Additional variance when rejecting
     19    ppStackOptions *options = args->data[1]; // Options
     20    pmConfig *config = args->data[2];   // Configuration
     21
     22    pmReadout *outRO = options->outRO;  // Output readout
     23    psVector *mask = options->inputMask; // Mask for inputs
     24    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
     25    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
    2326
    2427    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
    … …  
    3740    psArray *args = job->args;          // Arguments
    3841    ppStackThread *thread = args->data[0]; // Thread
    39     pmConfig *config = args->data[1];   // Configuration
    40     pmReadout *outRO = args->data[2];   // Output readout
    41     psVector *mask = args->data[3];     // Mask for inputs
    42     psArray *rejected = args->data[4];  // Rejected pixels
    43     psVector *weightings = args->data[5];  // Weighting (1/noise^2) for each image
    44     psVector *addVariance = args->data[6];  // Weighting (1/noise^2) for each image
     42    ppStackOptions *options = args->data[1]; // Options
     43    pmConfig *config = args->data[2];   // Configuration
     44
     45    pmReadout *outRO = options->outRO;  // Output readout
     46    psVector *mask = options->inputMask; // Mask for inputs
     47    psArray *rejected = options->rejected; // Rejected pixels
     48    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
     49    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
    4550
    4651    bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected,
    … …  
    227232        }
    228233
    229         pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]);
     234        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i],
     235                                             addVariance ? addVariance->data.F32[i] : NAN);
    230236        data->reject = psMemIncrRefCounter(rejected->data[i]);
    231237        stack->data[i] = data;

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

    r21260 r23352  
    1313#define FAKE_ROWS 4913
    1414
    15 float ppStackSourcesTransparency(const psArray *sourceLists, const pmFPAview *view, const pmConfig *config)
     15#ifdef TESTING
     16// Dump matches to a file
     17static void dumpMatches(const char *filename, // File to which to dump
     18                        int num,        // Number of inputs
     19                        psArray *matches, // Star matches
     20                        psVector *zp,   // Zero points
     21                        psVector *trans // Transparencies
     22                        )
     23{
     24    FILE *outMatches = fopen(filename, "w"); // Output matches
     25    psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star
     26    psVector *magErr = psVectorAlloc(num, PS_TYPE_F32); // Errors for each star
     27    for (int i = 0; i < matches->n; i++) {
     28        pmSourceMatch *match = matches->data[i]; // Match of interest
     29        psVectorInit(mag, NAN);
     30        psVectorInit(magErr, NAN);
     31        for (int j = 0; j < match->num; j++) {
     32            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     33                continue;
     34            }
     35            int index = match->image->data.U32[j]; // Image index
     36            mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index];
     37            if (trans) {
     38                mag->data.F32[index] -= trans->data.F32[index];
     39            }
     40            magErr->data.F32[index] = match->magErr->data.F32[j];
     41        }
     42        for (int j = 0; j < num; j++) {
     43            fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]);
     44        }
     45        fprintf(outMatches, "\n");
     46    }
     47    psFree(mag);
     48    psFree(magErr);
     49    fclose(outMatches);
     50    return;
     51}
     52#endif
     53
     54
     55float ppStackSourcesTransparency(const psArray *sourceLists, psVector *inputMask,
     56                                 const pmFPAview *view, const pmConfig *config)
    1657{
    1758    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);
    1862    PS_ASSERT_PTR_NON_NULL(view, NAN);
    1963    PS_ASSERT_PTR_NON_NULL(config, NAN);
    2064
    21 #ifdef TESTING
     65#if defined(TESTING) && 0
    2266    {
    2367        // Deliberately induce a major transparency difference
    … …  
    3781
    3882    float radius = psMetadataLookupF32(NULL, recipe, "ZP.RADIUS"); // Radius (pixels) for matching sources
    39     int iter = psMetadataLookupS32(NULL, recipe, "ZP.ITER"); // Maximum iterations
     83    int iter1 = psMetadataLookupS32(NULL, recipe, "ZP.ITER.1"); // Maximum iterations for pass 1
     84    int iter2 = psMetadataLookupS32(NULL, recipe, "ZP.ITER.2"); // Maximum iterations for pass 2
    4085    float tol = psMetadataLookupF32(NULL, recipe, "ZP.TOL"); // Tolerance for zero point iterations
    4186    int transIter = psMetadataLookupS32(NULL, recipe, "ZP.TRANS.ITER"); // Iterations for transparency
    4287    float transRej = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.REJ");// Rejection threshold for transparency
    4388    float transThresh = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.THRESH"); // Threshold for transparency
    44     float starRej = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ"); // Rejection threshold for stars
     89
     90    float starRej1 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ.1"); // Rejection threshold for stars
     91    float starSys1 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS.1"); // Estimated systematic error
     92    float starRej2 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ.2"); // Rejection threshold for stars
     93    float starSys2 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS.2"); // Estimated systematic error
     94
    4595    float starLimit = psMetadataLookupF32(NULL, recipe, "ZP.STAR.LIMIT"); // Limit on star rejection fraction
    46     float starSys = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS"); // Estimated systematic error
     96
     97    float fracMatch = psMetadataLookupF32(NULL, recipe, "ZP.MATCH"); // Fraction of images to match for star
    4798
    4899    psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms
    … …  
    63114        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
    64115
    65 #ifdef TESTING
     116#if defined(TESTING) && 0
    66117        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
    67118        pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
    68         pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sourceLists->data[i], NULL, NULL, psf, 5, 0, false, true);
     119        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sourceLists->data[i],
     120                                 NULL, NULL, psf, 5, 0, false, true);
    69121        psString name = NULL;
    70122        psStringAppend(&name, "start_%03d.fits", i);
     123        pmStackVisualPlotTestImage(fake->image, name);
    71124        psFits *fits = psFitsOpen(name, "w");
    72125        psFree(name);
    … …  
    108161    }
    109162
    110     psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
     163    psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
    111164    if (!matches) {
    112165        psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
    … …  
    114167        return NAN;
    115168    }
    116     psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
    117                                            transThresh, starRej, starSys); // Transparencies for each image
    118 
    119 #ifdef TESTING
    120     {
    121         // Dump the corrected magnitudes
    122         FILE *outMatches = fopen("source_match.dat", "w"); // Output matches
    123         psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star
    124         psVector *magErr = psVectorAlloc(num, PS_TYPE_F32); // Errors for each star
    125         for (int i = 0; i < matches->n; i++) {
    126             pmSourceMatch *match = matches->data[i]; // Match of interest
    127             psVectorInit(mag, NAN);
    128             psVectorInit(magErr, NAN);
    129             for (int j = 0; j < match->num; j++) {
    130                 if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
    131                     continue;
    132                 }
    133                 int index = match->image->data.U32[j]; // Image index
    134                 mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index] - trans->data.F32[index];
    135                 magErr->data.F32[index] = match->magErr->data.F32[j];
    136             }
    137             for (int j = 0; j < num; j++) {
    138                 fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]);
    139             }
    140             fprintf(outMatches, "\n");
    141         }
    142         psFree(mag);
    143         psFree(magErr);
    144         fclose(outMatches);
    145     }
    146 #endif
     169
     170#ifdef TESTING
     171    dumpMatches("source_match.dat", num, matches, zp, NULL);
     172#endif
     173
     174    psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
     175                                           iter2, starRej2, starSys2, starLimit,
     176                                           transIter, transRej, transThresh); // Transparencies for each image
     177    if (!trans) {
     178        psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
     179        return NAN;
     180    }
     181
     182#ifdef TESTING
     183    dumpMatches("source_mags.dat", num, matches, zp, trans);
     184#endif
     185
     186    for (int i = 0; i < trans->n; i++) {
     187        if (!isfinite(trans->data.F32[i])) {
     188            inputMask->data.U8[i] = PPSTACK_MASK_CAL;
     189        }
     190    }
    147191
    148192    // Save best matches SOMEWHERE for future photometry
    149193    // XXX this is a really poor output location; clean up the pmFPAfiles used in ppStack
    150194    pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT");
    151     psArray *sourcesBest = psArrayAllocEmpty (100);
    152 
    153     // XXX something of a hack: require at 2 detections or 1/2 of the max possible
    154     int minMatches = PS_MAX (2, 0.5*num);
     195    psArray *sourcesBest = psArrayAllocEmpty(matches->n);
     196
     197    // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible
     198    int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required
    155199    for (int i = 0; i < matches->n; i++) {
    156         pmSourceMatch *match = matches->data[i]; // Match of interest
    157         if (match->num < minMatches) continue;
    158 
    159         // We need to grab a single instance of this source: just take the first available
    160         int nImage = match->image->data.S32[0];
    161         int nIndex = match->index->data.S32[0];
    162         psArray *sources = sourceLists->data[nImage];
    163         pmSource *source = sources->data[nIndex];
    164        
    165         // stick this sample source on sourcesBest
    166         psArrayAdd (sourcesBest, 100, source);
    167     }
    168     psMetadataAdd (sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, "psphot sources", sourcesBest);
    169     psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis\n", sourcesBest->n);
    170     psFree (sourcesBest);
     200        pmSourceMatch *match = matches->data[i]; // Match of interest
     201        if (match->num < minMatches) {
     202            continue;
     203        }
     204
     205        // We need to grab a single instance of this source: just take the first available
     206        int image = match->image->data.S32[0]; // Index of image
     207        int index = match->index->data.S32[0]; // Index of source within image
     208        psArray *sources = sourceLists->data[image]; // Sources for image
     209        pmSource *source = sources->data[index]; // Source of interest
     210
     211        psArrayAdd(sourcesBest, sourcesBest->n, source);
     212    }
     213    psMetadataAdd(sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE,
     214                  "psphot sources", sourcesBest);
     215    psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
     216    psFree(sourcesBest);
    171217
    172218    psFree(matches);
    173     if (!trans) {
    174         psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
    175         return NAN;
    176     }
    177219
    178220    // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
    … …  
    182224    // We don't need to know the magnitude zero point for the filter, since it cancels out
    183225    for (int i = 0; i < num; i++) {
     226        if (!isfinite(trans->data.F32[i])) {
     227            continue;
     228        }
    184229        psArray *sources = sourceLists->data[i]; // Sources of interest
    185230        float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];

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

    r21477 r23352  
    99
    1010#include "ppStack.h"
     11#include "ppStackOptions.h"
     12#include "ppStackThread.h"
     13
    1114
    1215#define THREAD_WAIT 10000               // Microseconds to wait if thread is not available
    … …  
    3639    psFree(stack->threads);
    3740    for (int i = 0; i < stack->imageFits->n; i++) {
    38         psFitsClose(stack->imageFits->data[i]);
    39         psFitsClose(stack->maskFits->data[i]);
    40         psFitsClose(stack->varianceFits->data[i]);
     41        if (stack->imageFits->data[i]) {
     42            psFitsClose(stack->imageFits->data[i]);
     43        }
     44        if (stack->maskFits->data[i]) {
     45            psFitsClose(stack->maskFits->data[i]);
     46        }
     47        if (stack->varianceFits->data[i]) {
     48            psFitsClose(stack->varianceFits->data[i]);
     49        }
    4150        stack->imageFits->data[i] = stack->maskFits->data[i] = stack->varianceFits->data[i] = NULL;
    4251    }
    … …  
    4756}
    4857
    49 ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, const psArray *imageNames,
    50                                           const psArray *maskNames, const psArray *varianceNames,
    51                                           const psArray *covariances, const pmConfig *config)
    52 {
     58ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config)
     59{
     60    psAssert(options, "Require options");
     61    psAssert(config, "Require configuration");
     62
     63    const psArray *cells = options->cells; // Array of input cells
     64    const psArray *imageNames = options->imageNames; // Names of images to read
     65    const psArray *maskNames = options->maskNames; // Names of masks to read
     66    const psArray *varianceNames = options->varianceNames; // Names of variance maps to read
     67    const psArray *covariances = options->covariances; // Covariance matrices (already read)
     68
    5369    PS_ASSERT_ARRAY_NON_NULL(cells, NULL);
    5470    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL);
    … …  
    239255
    240256    {
    241         psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 6);
     257        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 3);
    242258        task->function = &ppStackReadoutInitialThread;
    243259        psThreadTaskAdd(task);
    … …  
    253269
    254270    {
    255         psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7);
     271        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 3);
    256272        task->function = &ppStackReadoutFinalThread;
    257273        psThreadTaskAdd(task);

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

    r18419 r23352  
    77#include <psmodules.h>
    88#include <ppStats.h>
     9#include <psphot.h>
    910
    1011#include "ppStack.h"
    1112
    12 static const char *cvsTag = "$Name: not supported by cvs2svn $";// CVS tag name
     13#ifndef PPSTACK_VERSION
     14#error "PPSTACK_VERSION is not set"
     15#endif
     16#ifndef PPSTACK_BRANCH
     17#error "PPSTACK_BRANCH is not set"
     18#endif
     19#ifndef PPSTACK_SOURCE
     20#error "PPSTACK_SOURCE is not set"
     21#endif
     22
     23#define xstr(s) str(s)
     24#define str(s) #s
    1325
    1426psString ppStackVersion(void)
    1527{
    16     psString version = NULL;            // Version, to return
    17     psStringAppend(&version, "%s-%s",PACKAGE_NAME,PACKAGE_VERSION);
    18     return version;
     28    char *value = NULL;
     29    psStringAppend(&value, "%s@%s", xstr(PPSTACK_BRANCH), xstr(PPSTACK_VERSION));
     30    return value;
     31}
     32
     33psString ppStackSource(void)
     34{
     35    return psStringCopy(xstr(PPSTACK_SOURCE));
    1936}
    2037
    2138psString ppStackVersionLong(void)
    2239{
    23     psString version = ppStackVersion(); // Version, to return
    24     psString tag = psStringStripCVS(cvsTag, "Name"); // CVS tag
    25     psStringAppend(&version, " (cvs tag %s) %s, %s", tag, __DATE__, __TIME__);
    26     psFree(tag);
     40    psString version = ppStackVersion();  // Version, to return
     41    psString source = ppStackSource();    // Source
     42
     43    psStringPrepend(&version, "ppStack ");
     44    psStringAppend(&version, " from %s, built %s, %s", source, __DATE__, __TIME__);
     45    psFree(source);
     46
     47#ifdef __OPTIMIZE__
     48    psStringAppend(&version, " optimised");
     49#else
     50    psStringAppend(&version, " unoptimised");
     51#endif
     52
    2753    return version;
    28 }
     54};
    2955
    3056
    31 void ppStackVersionMetadata(psMetadata *metadata)
     57bool ppStackVersionHeader(psMetadata *header)
    3258{
    33     PS_ASSERT_METADATA_NON_NULL(metadata,);
    34 
    35     psString pslib = psLibVersionLong();// psLib version
    36     psString psmodules = psModulesVersionLong(); // psModules version
    37     psString ppStats = ppStatsVersionLong(); // ppStats version
    38     psString ppStack = ppStackVersionLong(); // ppStack version
     59    PS_ASSERT_METADATA_NON_NULL(header, false);
    3960
    4061    psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now
    4162    psString timeString = psTimeToISO(time); // The time in an ISO string
    4263    psFree(time);
    43     psString head = NULL;               // Head string
    44     psStringAppend(&head, "ppStack processing at %s. Component information:", timeString);
     64    psString history = NULL;               // History string
     65    psStringAppend(&history, "ppStack at %s", timeString);
     66    psFree(timeString);
     67    psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history);
     68    psFree(history);
     69
     70    psLibVersionHeader(header);
     71    psModulesVersionHeader(header);
     72    psphotVersionHeader(header);
     73    ppStatsVersionHeader(header);
     74
     75    psString version = ppStackVersion(); // Software version
     76    psString source  = ppStackSource();  // Software source
     77
     78    psStringPrepend(&version, "ppStack version: ");
     79    psStringPrepend(&source, "ppStack source: ");
     80
     81    psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, version);
     82    psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, source);
     83
     84    psFree(version);
     85    psFree(source);
     86
     87    return true;
     88}
     89
     90void ppStackVersionPrint(void)
     91{
     92    psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now
     93    psString timeString = psTimeToISO(time); // The time in an ISO string
     94    psFree(time);
     95    psLogMsg("ppStack", PS_LOG_INFO, "ppStack at %s", timeString);
    4596    psFree(timeString);
    4697
    47     psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, "", head);
    48     psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, "", pslib);
    49     psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, "", psmodules);
    50     psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, "", ppStats);
    51     psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, "", ppStack);
     98    psString pslib = psLibVersionLong();// psLib version
     99    psString psmodules = psModulesVersionLong(); // psModules version
     100    psString psphot = psphotVersionLong(); // psphot version
     101    psString ppStats = ppStatsVersionLong(); // psastro version
     102    psString ppStack = ppStackVersionLong(); // ppStack version
    52103
    53     psFree(head);
     104    psLogMsg("ppStack", PS_LOG_INFO, "%s", pslib);
     105    psLogMsg("ppStack", PS_LOG_INFO, "%s", psmodules);
     106    psLogMsg("ppStack", PS_LOG_INFO, "%s", psphot);
     107    psLogMsg("ppStack", PS_LOG_INFO, "%s", ppStats);
     108    psLogMsg("ppStack", PS_LOG_INFO, "%s", ppStack);
     109
    54110    psFree(pslib);
    55111    psFree(psmodules);
     112    psFree(psphot);
    56113    psFree(ppStats);
    57114    psFree(ppStack);
Note: See TracChangeset for help on using the changeset viewer.