IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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.