IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 20995


Ignore:
Timestamp:
Dec 15, 2008, 1:18:30 PM (17 years ago)
Author:
Paul Price
Message:

Using new pmSourceMatch and pmSourceMatchRelphot to do the source matching and relative photometry. This should be more accurate than measuring the magnitude difference image by image.

Location:
trunk/ppStack/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStack.h

    r20753 r20995  
    152152
    153153
    154 /// Add sources to source lists, removing duplicates and solving for magnitude differences
     154/// Calculate transparency differences between images
    155155///
    156 /// Corrects the sources to have consistent magnitudes.  Returns the source lists.
    157 psArray *ppStackSourceListAdd(psArray *lists, // List to which to add, or NULL
    158                               psArray *sources, // Sources to add
    159                               const pmConfig *config // Configuration
    160     );
    161 
    162 /// Combine source lists to yield a set of unique sources.
    163 ///
    164 /// Corrects the sources to have consistent magnitudes where possible.  Returns the sources
    165 psArray *ppStackSourceListCombine(psArray *lists, // Source lists
    166                                   const pmConfig *config // Configuration
    167     );
    168 
    169 /// Print sources into a ds9 regions file
    170 void ppStackSourcesPrint(const psArray *sources // Sources to print
     156/// Corrects the source PSF photometry to a common system
     157bool ppStackSourcesTransparency(const psArray *sourceLists, // Sources for each input
     158                                const pmFPAview *view, // View to readout
     159                                const pmConfig *config // Configuration
    171160    );
    172161

  • trunk/ppStack/src/ppStackArguments.c

    r20497 r20995  
    164164    psMetadataAddBool(arguments, PS_LIST_TAIL, "-safe", 0,
    165165                      "Play safe with small numbers of pixels to combine?", false);
    166     psMetadataAddF32(arguments, PS_LIST_TAIL, "-source-radius", 0, "Source exclusion radius", NAN);
    167     psMetadataAddS32(arguments, PS_LIST_TAIL, "-source-iter", 0, "Source clipping iterations", 0);
    168     psMetadataAddF32(arguments, PS_LIST_TAIL, "-source-rej", 0, "Source clipping rejection level", NAN);
    169     psMetadataAddS32(arguments, PS_LIST_TAIL, "-source-min", 0, "Source minimum overlap", 0);
     166    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-radius", 0, "Radius (pixels) for matching sources", NAN);
     167    psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-iter", 0, "Maximum iterations for zero point", 0);
     168    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-tol", 0, "Tolerance for zero point iterations", NAN);
     169    psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-trans-iter", 0, "Iterations for transparency determination", 0);
     170    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-trans-rej", 0, "Rejection threshold for transparency determination", NAN);
     171    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-trans-thresh", 0, "Threshold for transparency determination", NAN);
     172    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-rej", 0, "Rejection threshold for stars", NAN);
     173    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-limit", 0, "Limit on star rejection fraction for successful iteration", NAN);
     174    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys", 0, "Estimated systematic error", NAN);
    170175    psMetadataAddBool(arguments, PS_LIST_TAIL, "-renorm", 0, "Renormalise variance maps?", false);
    171176    psMetadataAddStr(arguments, PS_LIST_TAIL, "-renorm-mean", 0,
    … …  
    238243    valueArgRecipeStr(arguments, recipe, "-mask-poor", "MASK.POOR", recipe);
    239244
    240     VALUE_ARG_RECIPE_FLOAT("-source-radius", "SOURCE.RADIUS", F32);
    241     VALUE_ARG_RECIPE_INT("-source-iter",     "SOURCE.ITER",   S32, 0);
    242     VALUE_ARG_RECIPE_FLOAT("-source-rej",    "SOURCE.REJ",    F32);
    243     VALUE_ARG_RECIPE_INT("-source-min",      "SOURCE.MIN",    S32, 0);
     245    VALUE_ARG_RECIPE_FLOAT("-zp-radius",       "ZP.RADIUS",       F32);
     246    VALUE_ARG_RECIPE_INT("-zp-iter",           "ZP.ITER",         S32, 0);
     247    VALUE_ARG_RECIPE_FLOAT("-zp-tol",          "ZP.TOL",          F32);
     248    VALUE_ARG_RECIPE_INT("-zp-trans-iter",     "ZP.TRANS.ITER",   S32, 0);
     249    VALUE_ARG_RECIPE_FLOAT("-zp-trans-rej",    "ZP.TRANS.REJ",    F32);
     250    VALUE_ARG_RECIPE_FLOAT("-zp-trans-thresh", "ZP.TRANS.THRESH", F32);
     251    VALUE_ARG_RECIPE_FLOAT("-zp-star-rej",     "ZP.STAR.REJ",     F32);
     252    VALUE_ARG_RECIPE_FLOAT("-zp-star-limit",   "ZP.STAR.LIMIT",   F32);
     253    VALUE_ARG_RECIPE_FLOAT("-zp-star-sys",     "ZP.STAR.SYS",     F32);
    244254
    245255    VALUE_ARG_RECIPE_INT("-psf-instances", "PSF.INSTANCES", S32, 0);

  • trunk/ppStack/src/ppStackCamera.c

    r20710 r20995  
    7373    bool haveWeights = false;           // Do we have weight maps?
    7474    bool havePSFs = false;              // Do we have PSFs?
    75     bool haveSources = (bool)psMetadataLookup(config->arguments, "PPSTACK.SOURCES"); // Have global sources?
    7675
    7776    psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info
    … …  
    195194
    196195        // Add the sources file
    197         if (!haveSources) {
    198             if (!sources || strlen(sources) == 0) {
    199                 if (havePSFs) {
    200                     psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SOURCES %d", i);
    201                     return false;
    202                 }
    203             } else {
    204                 if (!havePSFs) {
    205                     psWarning("SOURCES provided without PSF --- ignoring.");
    206                 } else {
    207                     psArray *sourcesFiles = psArrayAlloc(1); // Array of filenames for this FPA
    208                     sourcesFiles->data[0] = psMemIncrRefCounter(sources);
    209                     psMetadataAddArray(config->arguments, PS_LIST_TAIL, "SOURCES.FILENAMES", PS_META_REPLACE,
    210                                        "Filenames of SOURCES files", sourcesFiles);
    211                     psFree(sourcesFiles);
    212 
    213                     bool status;
    214                     pmFPAfile *sourcesFile = pmFPAfileBindFromArgs(&status, imageFile, config,
    215                                                                    "PPSTACK.INPUT.SOURCES",
    216                                                                    "SOURCES.FILENAMES");
    217                     if (!status) {
    218                         psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)",
    219                                 i, sources);
    220                         return false;
    221                     }
    222                     if (sourcesFile->type != PM_FPA_FILE_CMF) {
    223                         psError(PS_ERR_IO, true, "PPSTACK.INPUT.SOURCES is not of type CMF");
    224                         return false;
    225                     }
    226                 }
     196        {
     197            psArray *sourcesFiles = psArrayAlloc(1); // Array of filenames for this FPA
     198            sourcesFiles->data[0] = psMemIncrRefCounter(sources);
     199            psMetadataAddArray(config->arguments, PS_LIST_TAIL, "SOURCES.FILENAMES", PS_META_REPLACE,
     200                               "Filenames of SOURCES files", sourcesFiles);
     201            psFree(sourcesFiles);
     202
     203            bool status;
     204            pmFPAfile *sourcesFile = pmFPAfileBindFromArgs(&status, imageFile, config,
     205                                                           "PPSTACK.INPUT.SOURCES", "SOURCES.FILENAMES");
     206            if (!status) {
     207                psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)",
     208                        i, sources);
     209                return false;
     210            }
     211            if (sourcesFile->type != PM_FPA_FILE_CMF) {
     212                psError(PS_ERR_IO, true, "PPSTACK.INPUT.SOURCES is not of type CMF");
     213                return false;
    227214            }
    228215        }
    … …  
    269256    }
    270257
    271     if (haveSources) {
    272         // Global list of sources for use as stamps
    273         bool status = false;            // Found the file?
    274         if (havePSFs) {
    275             pmFPAfile *sources = pmFPAfileDefineFromArgs(&status, config, "PPSTACK.INPUT.SOURCES",
    276                                                          "PPSTACK.SOURCES");
    277             if (!status) {
    278                 psError(PS_ERR_IO, false, "Failed to load file definition PPSTACK.INPUT.SOURCES");
    279                 return false;
    280             }
    281             if (sources && sources->type != PM_FPA_FILE_CMF) {
    282                 psError(PS_ERR_IO, true, "PPSTACK.SOURCES is not of type CMF");
    283                 return false;
    284             }
    285         }
    286     }
    287 
    288258    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.NUM", 0, "Number of input files", i);
    289259    psMetadataAddBool(config->arguments, PS_LIST_TAIL, "HAVE.PSF", 0, "Have PSFs available?", havePSFs);
    290     psMetadataAddBool(config->arguments, PS_LIST_TAIL, "HAVE.SOURCES", 0, "Have global sources?",
    291                       haveSources);
    292260
    293261    // Output image

  • trunk/ppStack/src/ppStackLoop.c

    r20837 r20995  
    252252    // Preparation iteration: Load the sources, and get a target PSF model
    253253    psTrace("ppStack", 1, "Determining target PSF....\n");
    254     psArray *globalSources = NULL;      // Global list of sources for matching (haveSources = TRUE)
    255     psArray *indSources = psArrayAlloc(num); // Individual lists of sources for matching (haveSources = FALSE)
     254    psArray *sourceLists = psArrayAlloc(num); // Individual lists of sources for matching
    256255    pmPSF *targetPSF = NULL;            // Target PSF
    257     bool haveSources = psMetadataLookupBool(NULL, config->arguments, "HAVE.SOURCES"); // Global sources?
    258256    if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
    259257        pmFPAfileActivate(config->files, false, NULL);
    … …  
    269267        psMetadataItem *fileItem; // Item from iteration
    270268        psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
    271         psArray *sourceLists = NULL;    // Source lists for merging sources from multiple readouts
    272269        int numCols = 0, numRows = 0;   // Size of image
    273270        int index = 0;                  // Index for file
    … …  
    280277                psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
    281278                psFree(view);
    282                 psFree(globalSources);
    283                 psFree(indSources);
    284279                psFree(sourceLists);
    285280                psFree(fileIter);
    … …  
    298293                psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF.");
    299294                psFree(view);
    300                 psFree(globalSources);
    301                 psFree(indSources);
    302295                psFree(sourceLists);
    303296                psFree(fileIter);
    … …  
    310303            }
    311304
    312             if (!haveSources) {
    313                 pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
    314                 psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
    315                 if (!sources) {
    316                     psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
    317                     return NULL;
    318                 }
    319                 indSources->data[index] = psMemIncrRefCounter(sources);
    320 
    321 
    322                 // Calculate zero points if we don't trust the source lists
    323                 if (psMetadataLookupBool(NULL, recipe, "ZP")) {
    324                     psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
    325                     pmFPAfileActivate(config->files, false, NULL);
    326                     fileActivationSingle(config, convolveFiles, true, index);
    327                     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT",
    328                                                             index);// File of interest
    329                     pmFPAview *view = filesIterateDown(config);
    330                     if (!view) {
    331                         psFree(globalSources);
    332                         psFree(indSources);
    333                         psFree(targetPSF);
    334                         return false;
    335                     }
    336 
    337                     pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
    338 
    339                     if (!ppStackInputPhotometry(ro, sources, config)) {
    340                         psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
    341                         psFree(globalSources);
    342                         psFree(indSources);
    343                         psFree(targetPSF);
    344                         return false;
    345                     }
    346 
    347                     psFree(view);
    348                     if (!filesIterateUp(config)) {
    349                         psFree(globalSources);
    350                         psFree(indSources);
    351                         psFree(targetPSF);
    352                         return false;
    353                     }
    354                     pmFPAfileActivate(config->files, false, NULL);
    355                     fileActivation(config, prepareFiles, true);
    356                 }
    357 
    358 
    359 #ifdef TESTING
    360                 ppStackSourcesPrint(sources);
    361 #endif
    362                 sourceLists = ppStackSourceListAdd(sourceLists, sources, config);
    363                 if (!sourceLists) {
    364                     psError(PS_ERR_UNKNOWN, false, "Unable to add sources to list.");
    365                     psFree(view);
    366                     psFree(globalSources);
    367                     psFree(indSources);
    368                     psFree(fileIter);
    369                     psFree(psfs);
     305            pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
     306            psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
     307            if (!sources) {
     308                psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
     309                return NULL;
     310            }
     311            sourceLists->data[index] = psMemIncrRefCounter(sources);
     312
     313            // Re-do photometry if we don't trust the source lists
     314            if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
     315                psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
     316                pmFPAfileActivate(config->files, false, NULL);
     317                fileActivationSingle(config, convolveFiles, true, index);
     318                pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
     319                pmFPAview *view = filesIterateDown(config);
     320                if (!view) {
     321                    psFree(sourceLists);
     322                    psFree(targetPSF);
    370323                    return false;
    371324                }
     325
     326                pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
     327
     328                if (!ppStackInputPhotometry(ro, sources, config)) {
     329                    psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
     330                    psFree(sourceLists);
     331                    psFree(targetPSF);
     332                    return false;
     333                }
     334
     335                psFree(view);
     336                if (!filesIterateUp(config)) {
     337                    psFree(sourceLists);
     338                    psFree(targetPSF);
     339                    return false;
     340                }
     341                pmFPAfileActivate(config->files, false, NULL);
     342                fileActivation(config, prepareFiles, true);
    372343            }
    373344
    … …  
    375346        }
    376347        psFree(fileIter);
     348
     349        // Zero point calibration
     350        if (!ppStackSourcesTransparency(sourceLists, view, config)) {
     351            psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
     352            psFree(sourceLists);
     353            psFree(targetPSF);
     354            return false;
     355        }
     356
     357        exit(1);
    377358
    378359        targetPSF = ppStackPSF(config, numCols, numRows, psfs);
    … …  
    380361        if (!targetPSF) {
    381362            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
    382             psFree(globalSources);
    383             psFree(indSources);
     363            psFree(sourceLists);
    384364            psFree(view);
    385365            return false;
    … …  
    390370                         "Target PSF", targetPSF);
    391371        outChip->data_exists = true;
    392 
    393         if (haveSources) {
    394             // We want to hang on to the 'sources' even when its host FPA is blown away
    395             pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSTACK.INPUT.SOURCES");
    396             globalSources = psMetadataLookupPtr(NULL, sourcesRO->analysis, "PSPHOT.SOURCES");
    397             psMemIncrRefCounter(globalSources);
    398             if (!globalSources) {
    399                 psError(PS_ERR_UNKNOWN, true, "Unable to find sources.");
    400                 psFree(view);
    401                 return false;
    402             }
    403         } else {
    404             globalSources = ppStackSourceListCombine(sourceLists, config);
    405             psFree(sourceLists);
    406             if (!globalSources) {
    407                 psError(PS_ERR_UNKNOWN, false, "Unable to add sources to list.");
    408                 psFree(view);
    409                 psFree(psfs);
    410                 return false;
    411             }
    412 #ifdef TESTING
    413             ppStackSourcesPrint(globalSources);
    414 #endif
    415         }
    416372
    417373        psFree(view);
    … …  
    460416        pmFPAview *view = filesIterateDown(config);
    461417        if (!view) {
    462             psFree(globalSources);
    463             psFree(indSources);
     418            psFree(sourceLists);
    464419            psFree(targetPSF);
    465420            psFree(rng);
    … …  
    478433            psError(PS_ERR_UNKNOWN, true, "Sizes of input images don't match: %dx%d vs %dx%d",
    479434                    readout->image->numCols, readout->image->numRows, numCols, numRows);
    480             psFree(globalSources);
    481             psFree(indSources);
     435            psFree(sourceLists);
    482436            psFree(targetPSF);
    483437            psFree(rng);
    … …  
    489443        // Background subtraction, scaling and normalisation is performed automatically by the image matching
    490444        psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
    491         psArray *sources = haveSources ? globalSources : indSources->data[i]; // Sources for matching
    492445        psTimerStart("PPSTACK_MATCH");
    493446        if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i], &weightings->data.F32[i],
    494                           sources, targetPSF, rng, config)) {
     447                          sourceLists->data[i], targetPSF, rng, config)) {
    495448            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
    496449            inputMask->data.U8[i] = PPSTACK_MASK_MATCH;
    … …  
    532485        memDump("match");
    533486    }
    534     psFree(globalSources);
    535     psFree(indSources);
     487    psFree(sourceLists);
    536488    psFree(targetPSF);
    537489    psFree(rng);

  • trunk/ppStack/src/ppStackPhotometry.c

    r20777 r20995  
    2424    bool mdok;                          // Status of MD lookup
    2525
    26     float zpRadius = psMetadataLookupS32(&mdok, recipe, "ZP.RADIUS"); // Radius for ZP measurement
     26    float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement
    2727    if (!mdok) {
    28         psError(PS_ERR_UNKNOWN, true, "Unable to find ZP.RADIUS in recipe");
     28        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe");
    2929        return false;
    3030    }
    31     float zpSigma = psMetadataLookupF32(&mdok, recipe, "ZP.SIGMA"); // Gaussian sigma for photometry
     31    float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry
    3232    if (!mdok) {
    33         psError(PS_ERR_UNKNOWN, true, "Unable to find ZP.SIGMA in recipe");
     33        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe");
    3434        return false;
    3535    }
    36     float zpFrac = psMetadataLookupF32(&mdok, recipe, "ZP.FRAC"); // Fraction of good pixels for photometry
     36    float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry
    3737    if (!mdok) {
    38         psError(PS_ERR_UNKNOWN, true, "Unable to find ZP.FRAC in recipe");
     38        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe");
    3939        return false;
    4040    }

  • trunk/ppStack/src/ppStackSources.c

    r20944 r20995  
    11#include <stdio.h>
     2#include <math.h>
     3#include <string.h>
    24#include <pslib.h>
    35#include <psmodules.h>
    … …  
    57#include "ppStack.h"
    68
    7 #define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
    8                      PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
    9                      PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
    10 #define SOURCE_FAINTEST 50.0            // Faintest magnitude to consider
    11 #define SOURCES_MAX_LEAF 2              // Maximum number of points on a tree leaf
    12 
    13 //#define TESTING                         // Enable test output
    14 
    15 // Stack of sources
    16 typedef struct {
    17     psRegion *bounds;                   // Bounding box for sources
    18     psArray *sources;                   // Sources within region
    19     psVector *indices;                  // Indices from tree to sources
    20     psArray *duplicates;                // Duplicate sources (from overlaps)
    21     psTree *tree;                       // kd tree with sources
    22 } ppStackSourceList;
     9//#define TESTING                         // Enable debugging output
    2310
    2411
    25 // Extract coordinates from a source
    26 static inline void coordsFromSource(float *x, float *y, // Coordinates to return
    27                                     pmSource *source // Source of interest
    28     )
     12bool ppStackSourcesTransparency(const psArray *sourceLists, const pmFPAview *view, const pmConfig *config)
    2913{
    30     if (!source) {
    31         *x = NAN;
    32         *y = NAN;
    33     } else if (source->modelPSF) {
    34         *x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    35         *y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    36     } else {
    37         *x = source->peak->xf;
    38         *y = source->peak->yf;
    39     }
    40     return;
    41 }
     14    PS_ASSERT_ARRAY_NON_NULL(sourceLists, false);
     15    PS_ASSERT_PTR_NON_NULL(view, false);
     16    PS_ASSERT_PTR_NON_NULL(config, false);
    4217
    43 // Deallocator
    44 static void stackSourceListFree(ppStackSourceList *list)
    45 {
    46     psFree(list->bounds);
    47     psFree(list->sources);
    48     psFree(list->indices);
    49     psFree(list->duplicates);
    50     psFree(list->tree);
    51 }
    52 
    53 // Parse the sources into vectors for each coordinate, and a bounding box
    54 // Returns number of valid sources
    55 static long stackSourcesParse(psRegion **bounds, // Region to update with bounding box
    56                               psVector **x, psVector **y, // Coordinate vectors to return
    57                               psVector **indices, // Source indices to return
    58                               const psArray *sources // Input sources
    59     )
    60 {
    61     psAssert(bounds, "Must be given a region for bounding box");
    62     psAssert(x && y, "Must be given vectors");
    63     psAssert(indices, "Must be given indices");
    64     psAssert(sources, "Must be given sources");
    65 
    66     long numSources = sources->n;              // Number of sources
    67     *x = psVectorRecycle(*x, numSources, PS_TYPE_F32);
    68     *y = psVectorRecycle(*y, numSources, PS_TYPE_F32);
    69     *indices = psVectorRecycle(*indices, numSources, PS_TYPE_U32);
    70     float xMin = INFINITY, xMax = -INFINITY, yMin = INFINITY, yMax = -INFINITY; // Bounds of sources
    71     long num = 0;                       // Number of valid sources
    72     for (long i = 0; i < numSources; i++) {
    73         pmSource *source = sources->data[i]; // Source of interest
    74         if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) ||
    75             source->psfMag > SOURCE_FAINTEST) {
    76             continue;
    77         }
    78 
    79         float xSrc, ySrc;               // Coordinates of source
    80         coordsFromSource(&xSrc, &ySrc, source);
    81         if (xSrc < xMin) xMin = xSrc;
    82         if (xSrc > xMax) xMax = xSrc;
    83         if (ySrc < yMin) yMin = ySrc;
    84         if (ySrc > yMax) yMax = ySrc;
    85 
    86         (*x)->data.F32[num] = xSrc;
    87         (*y)->data.F32[num] = ySrc;
    88         (*indices)->data.U32[num] = i;
    89         num++;
    90     }
    91     (*x)->n = num;
    92     (*y)->n = num;
    93     (*indices)->n = num;
    94 
    95     if (*bounds) {
    96         (*bounds)->x0 = xMin;
    97         (*bounds)->x1 = xMax;
    98         (*bounds)->y0 = yMin;
    99         (*bounds)->y1 = yMax;
    100     } else {
    101         *bounds = psRegionAlloc(xMin, xMax, yMin, yMax);
    102     }
    103 
    104     psTrace("ppStack.sources", 8, "%ld sources: bounds are [%.2f:%.2f,%.2f:%.2f]\n",
    105             num, xMin, xMax, yMin, yMax);
    106 
    107     return num;
    108 }
    109 
    110 
    111 // Allocator for ppStackSourceList
    112 ppStackSourceList *ppStackSourceListAlloc(psArray *sources // Sources to add
    113     )
    114 {
    115     psAssert(sources, "Must be given sources");
    116 
    117     ppStackSourceList *list = psAlloc(sizeof(ppStackSourceList));
    118     psMemSetDeallocator(list, (psFreeFunc)stackSourceListFree);
    119 
    120     // Copy the sources onto a new array, so there's no confusion between the inputs and outputs
    121     long num = sources->n;              // Number of sources
    122     list->sources = psArrayAlloc(num);
    123     for (long i = 0; i < num; i++) {
    124         list->sources->data[i] = psMemIncrRefCounter(sources->data[i]);
    125     }
    126 
    127     list->bounds = NULL;
    128     list->indices = NULL;
    129     list->duplicates = psArrayAllocEmpty(0);
    130 
    131     psVector *x = NULL, *y = NULL;      // Source coordinates
    132     stackSourcesParse(&list->bounds, &x, &y, &list->indices, sources);
    133     list->tree = psTreePlant(2, SOURCES_MAX_LEAF, x, y); // kd Tree with sources
    134     psFree(x);
    135     psFree(y);
    136 
    137     return list;
    138 }
    139 
    140 // Extend a list of sources
    141 bool ppStackSourceListExtend(ppStackSourceList *list, // List to extend
    142                              const psArray *newSources, // Sources to add
    143                              const psArray *newDups // Duplicate sources to add
    144     )
    145 {
    146     psAssert(list, "Must be given a list");
    147     psAssert(newSources, "Must be given sources");
    148 
     18#ifdef TESTING
    14919    {
    150         psArray *oldSources = list->sources;// Old list of sources
    151         long numOld = oldSources->n;        // Number of old sources
    152         long numNew = newSources->n;        // Number of new sources
    153 
    154         list->sources = oldSources = psArrayRealloc(oldSources, numOld + numNew);
    155         for (long i = 0; i < numNew; i++) {
    156             psArrayAdd(oldSources, 1, newSources->data[i]);
    157         }
    158         psFree(list->tree);
    159 
    160         psVector *x = NULL, *y = NULL;      // Source coordinates
    161         stackSourcesParse(&list->bounds, &x, &y, &list->indices, list->sources);
    162         list->tree = psTreePlant(2, SOURCES_MAX_LEAF, x, y); // kd Tree with sources
    163         psFree(x);
    164         psFree(y);
    165     }
    166 
    167     if (newDups) {
    168         psArray *oldDups = list->duplicates; // Old list of duplicates
    169         long numOld = oldDups ? oldDups->n : 0; // Number of old duplicates
    170         long numNew = newDups->n;       // Number of new duplicates
    171 
    172         list->duplicates = oldDups = psArrayRealloc(oldDups, numOld + numNew);
    173         for (long i = 0; i < numNew; i++) {
    174             psArrayAdd(oldDups, 1, newDups->data[i]);
     20        // Deliberately induce a major transparency difference
     21        psArray *sources = sourceLists->data[1]; // Sources to correct
     22        for (int i = 0; i < sources->n; i++) {
     23            pmSource *source = sources->data[i]; // Source of interest
     24            if (!source) {
     25                continue;
     26            }
     27            source->psfMag += 1.0;
    17528        }
    17629    }
    177 
    178     return list;
    179 }
    180 
    181 // Compare an array of sources with a list
    182 // Return number inside the list
    183 static long stackSourceListCompare(ppStackSourceList **merge, // Merge target, modified
    184                                    psArray *lists, // Array of source lists
    185                                    int listIndex, // Index of list to compare with sources
    186                                    psArray **sourcesPtr, // Sources to compare with list, modified
    187                                    float radius, // Matching radius
    188                                    int minOverlap, // Minimum number of sources
    189                                    int iter, // Clipping iterations
    190                                    float rej // Clipping rejection level
    191     )
    192 {
    193     psAssert(merge, "Expected a source list pointer");
    194     psAssert(radius > 0, "Silly radius: %f", radius);
    195     psAssert(lists, "Expected an array of source lists");
    196     psAssert(sourcesPtr && *sourcesPtr, "Expected a source array");
    197 
    198     ppStackSourceList *list = lists->data[listIndex]; // List of interest
    199     if (!list) {
    200         return 0;
    201     }
    202 
    203     psArray *sources = *sourcesPtr;     // Sources to compare with list
    204     psVector *x = NULL, *y = NULL;      // Coordinates of sources
    205     psVector *indices = NULL;           // Indices for sources
    206     psRegion *bounds = NULL;            // Bounding box for sources
    207     long numSources = stackSourcesParse(&bounds, &x, &y, &indices, sources); // Number of (good) sources
    208 
    209     if (bounds->x0 > list->bounds->x1 || bounds->x1 < list->bounds->x0 ||
    210         bounds->y0 > list->bounds->y1 || bounds->y1 < list->bounds->y0) {
    211         // Bounds don't overlap, so the sources don't either
    212         psTrace("ppStack.sources", 7, "Bounds don't overlap\n");
    213         psFree(x);
    214         psFree(y);
    215         psFree(indices);
    216         psFree(bounds);
    217         return 0;
    218     }
    219     psFree(bounds);
    220 
    221     long numIn = 0;                     // Number of sources in list
    222     long numOut = 0;                    // Number of sources outside list
    223 
    224 #ifdef TESTING
    225     static int fileNum = 0;             // Number of file
    226     psString filename = NULL;   // Name of file
    227     psStringAppend(&filename, "source_match_%d.txt", fileNum++);
    228     FILE *file = fopen(filename, "w"); // File to write
    229     psFree(filename);
    230     fprintf(file, "# x1 y1 x2 y2 mag1 mag1err mag2 mag2err flag1 flag2\n");
    23130#endif
    232 
    233     psVector *magDiff = psVectorAlloc(numSources, PS_TYPE_F32); // Magnitude differences
    234     psVector *coords = psVectorAlloc(2, PS_TYPE_F32); // Coordinates of source
    235     psArray *outside = psArrayAllocEmpty(numSources); // Sources that don't correlate
    236     psArray *inside = psArrayAllocEmpty(numSources); // Sources that do correlate
    237     for (long i = 0; i < numSources; i++) {
    238         coords->data.F32[0] = x->data.F32[i];
    239         coords->data.F32[1] = y->data.F32[i];
    240         long match = psTreeNearestWithin(list->tree, coords, radius); // Index of match
    241         if (match < 0) {
    242             numOut++;
    243             psArrayAdd(outside, 1, sources->data[indices->data.U32[i]]);
    244             continue;
    245         }
    246 
    247         pmSource *listSource = list->sources->data[list->indices->data.U32[match]]; // Source in list
    248         pmSource *source = sources->data[indices->data.U32[i]]; // Source of interest
    249         float listMag = listSource->psfMag; // Magnitude of source in list
    250         float sourceMag = source->psfMag; // Magnitude of source
    251         if (!(listSource->mode & SOURCE_MASK) && !(source->mode & SOURCE_MASK) &&
    252             isfinite(listMag) && isfinite(sourceMag)) {
    253 #ifdef TESTING
    254             float xList, yList;         // Coordinates from source on list
    255             coordsFromSource(&xList, &yList, listSource);
    256             fprintf(file, "%f %f %f %f %f %f %f %f %.4x %.4x\n", x->data.F32[i], y->data.F32[i], xList, yList,
    257                     listSource->psfMag, listSource->errMag, source->psfMag, source->errMag,
    258                     listSource->mode, source->mode);
    259 #endif
    260             magDiff->data.F32[numIn] = listSource->psfMag - source->psfMag;
    261             psArrayAdd(inside, 1, source);
    262             numIn++;
    263         }
    264     }
    265     magDiff->n = numIn;
    266 
    267 #ifdef TESTING
    268     fclose(file);
    269 #endif
    270 
    271     psFree(coords);
    272     psFree(x);
    273     psFree(y);
    274     psFree(indices);
    275 
    276     psTrace("ppStack.sources", 7, "%ld sources (vs %d) overlap list %d\n", numIn, minOverlap, listIndex);
    277 
    278     if (numIn > minOverlap) {
    279         // There's sufficient overlap --- calculate the magnitude difference
    280         psStats *stats = psStatsAlloc(iter == 0 ? (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV) :
    281                                       (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV));
    282         stats->clipSigma = rej;
    283         stats->clipIter = iter;
    284         if (!psVectorStats(stats, magDiff, NULL, NULL, 0)) {
    285             psError(PS_ERR_UNKNOWN, false, "Unable to determine mean magnitude difference");
    286             psFree(magDiff);
    287             psFree(stats);
    288             return -1;
    289         }
    290         float meanDiff = iter == 0 ? stats->sampleMean : stats->clippedMean; // Mean magnitude difference
    291         float stdev = iter == 0 ? stats->sampleStdev : stats->clippedStdev; // Standard deviation
    292         long numStats = iter == 0 ? numIn : stats->clippedNvalues; // Number used for statistics
    293         psLogMsg("ppStack.sources", PS_LOG_INFO, "Magnitude difference from %ld overlaps is %f +/- %f\n",
    294                  numStats, meanDiff, stdev);
    295         psFree(stats);
    296 
    297         // Merge our sources into the list
    298 
    299         if (*merge) {
    300             // These sources act as a bridge from one list to another
    301             // Correct magnitudes to be on the system defined by the first source list
    302             // Correct the list to 'sources' (because 'sources' has been corrected to the merge target)
    303             long num = list->sources->n; // Number of sources
    304             for (long j = 0; j < num; j++) {
    305                 pmSource *source = list->sources->data[j];
    306                 source->psfMag -= meanDiff;
    307             }
    308             for (long j = 0; j < list->duplicates->n; j++) {
    309                 pmSource *source = list->duplicates->data[j];
    310                 source->psfMag -= meanDiff;
    311             }
    312 
    313             // Suck sources from this list into the first one we found.
    314             psTrace("ppStack.sources", 4, "Bridge: Merging lists.");
    315 
    316             ppStackSourceListExtend(*merge, list->sources, NULL);
    317             psFree(list);
    318             lists->data[listIndex] = NULL;
    319         } else {
    320             // Correct 'sources' magnitudes to match the list magnitudes
    321             for (long j = 0; j < numOut; j++) {
    322                 pmSource *source = outside->data[j]; // New source
    323                 source->psfMag += meanDiff;
    324             }
    325             for (long j = 0; j < numIn; j++) {
    326                 pmSource *source = inside->data[j]; // New source
    327                 source->psfMag += meanDiff;
    328             }
    329             // Put the sources that didn't correlate into the list
    330             psTrace("ppStack.sources", 4, "Overlap: Extending list");
    331             ppStackSourceListExtend(list, outside, inside);
    332             *merge = list;
    333         }
    334 
    335         // We only need to consider sources outside the list, since the lists are (mostly) disjoint
    336         psFree(*sourcesPtr);
    337         *sourcesPtr = outside;
    338     } else {
    339         psFree(outside);
    340     }
    341     psFree(inside);
    342     psFree(magDiff);
    343 
    344     return numIn;
    345 }
    346 
    347 
    348 psArray *ppStackSourceListAdd(psArray *lists, psArray *sources, const pmConfig *config)
    349 {
    350     PS_ASSERT_PTR_NON_NULL(sources, NULL);
    351     PS_ASSERT_PTR_NON_NULL(config, NULL);
    352 
    353     if (!lists) {
    354         lists = psArrayAllocEmpty(1);
    355     }
    356     if (lists->n == 0) {
    357         ppStackSourceList *list = ppStackSourceListAlloc(sources);
    358         psArrayAdd(lists, 1, list);
    359         psFree(list);                   // Drop reference
    360         return lists;
    361     }
    362 
    363     // Read recipe values
    364     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    365     psAssert(recipe, "We've thrown an error on this before.");
    366     float radius = psMetadataLookupF32(NULL, recipe, "SOURCE.RADIUS"); // Exclusion radius
    367     int minOverlap = psMetadataLookupS32(NULL, recipe, "SOURCE.MIN"); // Minimum number
    368     int iter = psMetadataLookupS32(NULL, recipe, "SOURCE.ITER"); // Clipping iterations
    369     float rej = psMetadataLookupF32(NULL, recipe, "SOURCE.REJ"); // Clipping rejection
    370 
    371     // Look for overlaps with existing lists
    372     int numLists = lists->n;            // Number of source lists on the stack
    373     ppStackSourceList *merge = NULL;    // Merged source list
    374 
    375     psMemIncrRefCounter(sources);       // So we can play with the pointer
    376     for (int i = 0; i < numLists; i++) {
    377         long numIn = stackSourceListCompare(&merge, lists, i, &sources, radius, minOverlap,
    378                                             iter, rej); // Number of overlapping sources
    379         if (numIn == -1) {
    380             psError(PS_ERR_UNKNOWN, false, "Unable to compare sources.");
    381             psFree(sources);
    382             psFree(lists);
    383             return NULL;
    384         }
    385     }
    386     psFree(sources);
    387 
    388     if (!merge) {
    389         // The source list is disjoint, so throw them into a new list
    390         psTrace("ppStack.sources", 4, "Disjoint source: Creating new list");
    391         ppStackSourceList *list = ppStackSourceListAlloc(sources);
    392         psArrayAdd(lists, 1, list);
    393         psFree(list);                   // Drop reference
    394     }
    395 
    396     return lists;
    397 }
    398 
    399 
    400 
    401 psArray *ppStackSourceListCombine(psArray *lists, const pmConfig *config)
    402 {
    403     PS_ASSERT_ARRAY_NON_NULL(lists, NULL);
    404     PS_ASSERT_PTR_NON_NULL(config, NULL);
    40531
    40632    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    40733    psAssert(recipe, "We've thrown an error on this before.");
    408     float radius = psMetadataLookupF32(NULL, recipe, "SOURCE.RADIUS"); // Exclusion radius
    40934
    410     int numLists = 0;                   // Number of lists
    411     ppStackSourceList *firstList = NULL; // Last list
    412     for (int i = 0; i < lists->n; i++) {
    413         ppStackSourceList *list = lists->data[i]; // List of interest
    414         if (list) {
    415             numLists++;
    416             if (!firstList) {
    417                 firstList = list;
    418             }
    419         }
     35    float radius = psMetadataLookupF32(NULL, recipe, "ZP.RADIUS"); // Radius (pixels) for matching sources
     36    int iter = psMetadataLookupS32(NULL, recipe, "ZP.ITER"); // Maximum iterations
     37    float tol = psMetadataLookupF32(NULL, recipe, "ZP.TOL"); // Tolerance for zero point iterations
     38    int transIter = psMetadataLookupS32(NULL, recipe, "ZP.TRANS.ITER"); // Iterations for transparency
     39    float transRej = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.REJ");// Rejection threshold for transparency
     40    float transThresh = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.THRESH"); // Threshold for transparency
     41    float starRej = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ"); // Rejection threshold for stars
     42    float starLimit = psMetadataLookupF32(NULL, recipe, "ZP.STAR.LIMIT"); // Limit on star rejection fraction
     43    float starSys = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS"); // Estimated systematic error
     44
     45    psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms
     46    if (!airmassZP) {
     47        psError(PS_ERR_UNKNOWN, false, "Unable to find ZP.AIRMASS in recipe.");
     48        return false;
    42049    }
    42150
    422     if (!firstList || numLists == 0) {
    423         psError(PS_ERR_UNKNOWN, true, "No lists found.");
    424         return NULL;
     51    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
     52    psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n);
     53
     54    psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Zero points for each image
     55    const char *filter = NULL;          // Filter name
     56    float airmassTerm = NAN;            // Airmass term
     57    float sumExpTime = 0.0;             // Sum of the exposure time
     58    for (int i = 0; i < num; i++) {
     59        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
     60        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
     61
     62        float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time
     63        float airmass = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.AIRMASS"); // Airmass
     64        const char *expFilter = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.FILTER"); // Filter name
     65        if (!isfinite(exptime) || exptime == 0 || !isfinite(airmass) || airmass == 0 ||
     66            !expFilter || strlen(expFilter) == 0) {
     67            psError(PS_ERR_UNEXPECTED_NULL, false,
     68                    "Unable to find exposure time (%f), airmass (%f) or filter (%s)",
     69                    exptime, airmass, expFilter);
     70            psFree(zp);
     71            return false;
     72        }
     73
     74        if (!filter) {
     75            filter = expFilter;
     76            airmassTerm = psMetadataLookupF32(NULL, airmassZP, filter);
     77            if (!isfinite(airmassTerm)) {
     78                psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     79                        "Unable to find airmass term (ZP.AIRMASS) for filter %s", filter);
     80                psFree(zp);
     81                return false;
     82            }
     83        } else if (strcmp(filter, expFilter) != 0) {
     84            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Filters don't match: %s vs %s", filter, expFilter);
     85            psFree(zp);
     86            return false;
     87        }
     88
     89        zp->data.F32[i] = airmassTerm * airmass + 2.5 * log10(exptime);
     90        sumExpTime += exptime;
    42591    }
    42692
    427 #if 0
    428     if (numLists == 1) {
    429         // Trivial case
    430         return psMemIncrRefCounter(firstList->sources);
     93    psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
     94    if (!matches) {
     95        psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
     96        psFree(zp);
     97        return false;
     98    }
     99    psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
     100                                           transThresh, starRej, starSys); // Transparencies for each image
     101
     102#ifdef TESTING
     103    {
     104        // Dump the corrected magnitudes
     105        FILE *outMatches = fopen("source_match.dat", "w"); // Output matches
     106        psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star
     107        for (int i = 0; i < matches->n; i++) {
     108            pmSourceMatch *match = matches->data[i]; // Match of interest
     109            psVectorInit(mag, NAN);
     110            for (int j = 0; j < match->num; j++) {
     111                if (match->mask->data.PS_TYPE_MASK_DATA[j]) {
     112                    continue;
     113                }
     114                int index = match->image->data.U32[j]; // Image index
     115                mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index] - trans->data.F32[index];
     116            }
     117            for (int j = 0; j < num; j++) {
     118                fprintf(outMatches, "%f ", mag->data.F32[j]);
     119            }
     120            fprintf(outMatches, "\n");
     121        }
     122        psFree(mag);
     123        fclose(outMatches);
    431124    }
    432125#endif
    433126
    434     for (int i = lists->n - 1; lists->data[i] != firstList; i--) {
    435         // There may be some small overlap (we know it's less than SOURCE.MIN), so attempt to correct
    436         // magnitudes of these first
    437         ppStackSourceList *lastList = lists->data[i]; // Last source list in the array
    438         if (!lastList) {
    439             continue;
    440         }
    441         psArray *sources = psMemIncrRefCounter(lastList->sources); // Sources of interest
    442         ppStackSourceList *merge = NULL; // Merge target
    443         for (int j = 0; j < i; j++) {
    444             ppStackSourceList *list = lists->data[i]; // List to compare
    445             if (!list) {
     127    psFree(matches);
     128    if (!trans) {
     129        psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
     130        return false;
     131    }
     132
     133    // M = m + c0 + c1 * airmass + 2.5log(t) + transparency
     134    // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
     135    // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
     136    // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
     137    // We don't need to know the magnitude zero point for the filter, since it cancels out
     138    for (int i = 0; i < num; i++) {
     139        psArray *sources = sourceLists->data[i]; // Sources of interest
     140        float magCorr = airmassTerm + 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];
     141        psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr);
     142
     143        for (int j = 0; j < sources->n; j++) {
     144            pmSource *source = sources->data[j]; // Source of interest
     145            if (!source) {
    446146                continue;
    447147            }
    448             // Note: setting iterations = 0, so no clipping (figure not enough values to do clipping)
    449             long numIn = stackSourceListCompare(&merge, lists, j, &sources, radius, 0,
    450                                                 1, NAN); // Number of overlapping sources
    451             if (numIn < 0) {
    452                 psError(PS_ERR_UNKNOWN, false, "Unable to compare sources.");
    453                 psFree(sources);
    454                 return NULL;
    455             }
    456             if (numIn == 0) {
    457                 // It doesn't match anything at all.  Merge it in to the first list.
    458                 psTrace("ppStack.sources", 4, "Disjoint list: Merging into first");
    459                 ppStackSourceListExtend(firstList, list->sources, NULL);
    460             }
     148            source->psfMag += magCorr;
    461149        }
    462         psFree(sources);
    463150    }
    464 
    465     return psMemIncrRefCounter(firstList->sources);
    466 }
     151    psFree(trans);
    467152
    468153
    469 
    470 void ppStackSourcesPrint(const psArray *sources)
    471 {
    472     static int num = 0;                 // Number of source array
    473     psString filename = NULL;           // Name of file
    474     psStringAppend(&filename, "sources_%d.reg", num++);
    475 
    476     const char *goodColour = "green";   // Colour for good sources
    477     const char *badColour = "red";      // Colour for bad sources
    478     float radius = 5;                   // Radius for circle
    479 
    480     FILE *file = fopen(filename, "w");  // File to which to write
    481     psFree(filename);
    482     for (int i = 0; i < sources->n; i++) {
    483         pmSource *source = sources->data[i]; // Source of interest
    484         if (!source) {
    485             continue;
     154#ifdef TESTING
     155    // Double check: all transparencies should be zero
     156    {
     157        psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
     158        if (!matches) {
     159            psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
     160            psFree(zp);
     161            return false;
    486162        }
    487 
    488         float x, y;                     // Source coordinates
    489         coordsFromSource(&x, &y, source);
    490 
    491         bool bad = (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) ||
    492             source->psfMag > SOURCE_FAINTEST; // Is this a bad source?
    493 
    494         fprintf(file, "image; circle(%f,%f,%f) # color = %s\n", x + 1.0, y + 1.0, radius,
    495                 bad ? badColour : goodColour);
     163        psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
     164                                               transThresh, starRej, starSys);
     165        psFree(trans);
    496166    }
    497     fclose(file);
    498 
    499     return;
     167#endif
     168    return true;
    500169}
Note: See TracChangeset for help on using the changeset viewer.