IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 23190 for trunk/ppStack


Ignore:
Timestamp:
Mar 5, 2009, 9:19:00 AM (17 years ago)
Author:
Paul Price
Message:

Bad images can be identified as those which have a NAN transparency correction. Extended the inputMask to cover these as well.

Location:
trunk/ppStack/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStack.h

    r23171 r23190  
    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;
    … …  
    161162/// Corrects the source PSF photometry to a common system.  Return the sum of the exposure times.
    162163float ppStackSourcesTransparency(const psArray *sourceLists, // Sources for each input
     164                                 psVector *inputMask, // Indicates bad input
    163165                                 const pmFPAview *view, // View to readout
    164166                                 const pmConfig *config // Configuration

  • trunk/ppStack/src/ppStackLoop.c

    r23171 r23190  
    257257    pmPSF *targetPSF = NULL;            // Target PSF
    258258    float sumExposure = NAN;            // Sum of exposure times
     259    psVector *inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs
     260    psVectorInit(inputMask, 0);
    259261    if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
    260262        pmFPAfileActivate(config->files, false, NULL);
    … …  
    283285                psFree(fileIter);
    284286                psFree(psfs);
     287                psFree(inputMask);
    285288                return false;
    286289            }
    … …  
    299302                psFree(fileIter);
    300303                psFree(psfs);
     304                psFree(inputMask);
    301305                return false;
    302306            }
    … …  
    324328                    psFree(sourceLists);
    325329                    psFree(targetPSF);
    326                     return false;
     330                    psFree(inputMask);
     331                   return false;
    327332                }
    328333
    … …  
    333338                    psFree(sourceLists);
    334339                    psFree(targetPSF);
     340                    psFree(inputMask);
    335341                    return false;
    336342                }
    … …  
    340346                    psFree(sourceLists);
    341347                    psFree(targetPSF);
     348                    psFree(inputMask);
    342349                    return false;
    343350                }
    … …  
    351358
    352359        // Zero point calibration
    353         sumExposure = ppStackSourcesTransparency(sourceLists, view, config);
     360        sumExposure = ppStackSourcesTransparency(sourceLists, inputMask, view, config);
    354361        if (!isfinite(sumExposure) || sumExposure <= 0) {
    355362            psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
    356363            psFree(sourceLists);
    357364            psFree(targetPSF);
     365            psFree(inputMask);
    358366            return false;
    359367        }
    … …  
    366374            psFree(sourceLists);
    367375            psFree(view);
     376            psFree(inputMask);
    368377            return false;
    369378        }
    … …  
    409418    int numGood = 0;                    // Number of good frames
    410419    int numCols = 0, numRows = 0;       // Size of image
    411     psVector *inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs
    412     psVectorInit(inputMask, 0);
    413420    psVector *matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching
    414421    psVectorInit(matchChi2, NAN);
    … …  
    419426    psArray *covariances = psArrayAlloc(num); // Covariance matrices
    420427    for (int i = 0; i < num; i++) {
     428        if (inputMask->data.U8[i]) {
     429            continue;
     430        }
    421431        psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num);
    422432        pmFPAfileActivate(config->files, false, NULL);

  • trunk/ppStack/src/ppStackSources.c

    r21260 r23190  
    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
    … …  
    63107        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
    64108
    65 #ifdef TESTING
     109#if defined(TESTING) && 0
    66110        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
    67111        pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
    … …  
    114158        return NAN;
    115159    }
     160
     161#ifdef TESTING
     162    dumpMatches("source_match.dat", num, matches, zp, NULL);
     163#endif
     164
    116165    psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
    117166                                           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
     167    if (!trans) {
     168        psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
     169        return NAN;
     170    }
     171
     172#ifdef TESTING
     173    dumpMatches("source_mags.dat", num, matches, zp, trans);
     174#endif
     175
     176    for (int i = 0; i < trans->n; i++) {
     177        if (!isfinite(trans->data.F32[i])) {
     178            inputMask->data.U8[i] = PPSTACK_MASK_CAL;
     179        }
     180    }
    147181
    148182    // Save best matches SOMEWHERE for future photometry
    … …  
    154188    int minMatches = PS_MAX (2, 0.5*num);
    155189    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);
     190        pmSourceMatch *match = matches->data[i]; // Match of interest
     191        if (match->num < minMatches) continue;
     192
     193        // We need to grab a single instance of this source: just take the first available
     194        int nImage = match->image->data.S32[0];
     195        int nIndex = match->index->data.S32[0];
     196        psArray *sources = sourceLists->data[nImage];
     197        pmSource *source = sources->data[nIndex];
     198
     199        // stick this sample source on sourcesBest
     200        psArrayAdd (sourcesBest, 100, source);
    167201    }
    168202    psMetadataAdd (sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, "psphot sources", sourcesBest);
    … …  
    171205
    172206    psFree(matches);
    173     if (!trans) {
    174         psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
    175         return NAN;
    176     }
    177207
    178208    // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
    … …  
    182212    // We don't need to know the magnitude zero point for the filter, since it cancels out
    183213    for (int i = 0; i < num; i++) {
     214        if (!isfinite(trans->data.F32[i])) {
     215            continue;
     216        }
    184217        psArray *sources = sourceLists->data[i]; // Sources of interest
    185218        float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];

  • trunk/ppStack/src/ppStackThread.c

    r21477 r23190  
    3636    psFree(stack->threads);
    3737    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]);
     38        if (stack->imageFits->data[i]) {
     39            psFitsClose(stack->imageFits->data[i]);
     40        }
     41        if (stack->maskFits->data[i]) {
     42            psFitsClose(stack->maskFits->data[i]);
     43        }
     44        if (stack->varianceFits->data[i]) {
     45            psFitsClose(stack->varianceFits->data[i]);
     46        }
    4147        stack->imageFits->data[i] = stack->maskFits->data[i] = stack->varianceFits->data[i] = NULL;
    4248    }
Note: See TracChangeset for help on using the changeset viewer.