IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21106


Ignore:
Timestamp:
Jan 12, 2009, 1:03:29 PM (17 years ago)
Author:
Paul Price
Message:

Adding more debugging output. Think I've solved the problem of the poor photometry: the background being added in to the fake images produces the bad photometry. Without it, the photometry s.d. goes from 0.08 mag to 0.02 mag. Want to replace it with the proper psphot background model (not done yet). These files are currently in an intermediate state.

Location:
branches/pap_branch_20090108/ppStack/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/pap_branch_20090108/ppStack/src/ppStackLoop.c

    r21101 r21106  
    348348        psFree(fileIter);
    349349
     350        // Generate target PSF
     351        targetPSF = ppStackPSF(config, numCols, numRows, psfs);
     352        psFree(psfs);
     353        if (!targetPSF) {
     354            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
     355            psFree(sourceLists);
     356            psFree(view);
     357            return false;
     358        }
     359        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN,
     360                         "Target PSF for stack", targetPSF);
     361
    350362        // Zero point calibration
    351363        sumExposure = ppStackSourcesTransparency(sourceLists, view, config);
    … …  
    354366            psFree(sourceLists);
    355367            psFree(targetPSF);
    356             return false;
    357         }
    358 
    359         targetPSF = ppStackPSF(config, numCols, numRows, psfs);
    360         psFree(psfs);
    361         if (!targetPSF) {
    362             psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
    363             psFree(sourceLists);
    364             psFree(view);
    365368            return false;
    366369        }
    … …  
    450453        psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
    451454        psTimerStart("PPSTACK_MATCH");
     455
     456
     457        // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
     458        // Passing common set of sources, rather than individual
     459
    452460        if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i], &weightings->data.F32[i],
    453                           sourceLists->data[i], targetPSF, rng, config)) {
     461                          sourceLists->data[0], targetPSF, rng, config)) {
    454462            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
    455463            inputMask->data.U8[i] = PPSTACK_MASK_MATCH;

  • branches/pap_branch_20090108/ppStack/src/ppStackMatch.c

    r21096 r21106  
    4747#endif
    4848
     49// Get coordinates from a source
     50static void coordsFromSource(float *x, float *y, const pmSource *source)
     51{
     52    assert(x && y);
     53    assert(source);
     54
     55    if (source->modelPSF) {
     56        *x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
     57        *y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
     58    } else {
     59        *x = source->peak->xf;
     60        *y = source->peak->yf;
     61    }
     62    return;
     63}
     64
    4965static psArray *stackSourcesFilter(psArray *sources, // Source list to filter
    5066                                   int exclusion // Exclusion zone, pixels
    … …  
    6480            continue;
    6581        }
    66         x->data.F32[numGood] = source->peak->xf;
    67         y->data.F32[numGood] = source->peak->yf;
     82        coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source);
    6883        numGood++;
    6984    }
    … …  
    8095            continue;
    8196        }
    82         coords->data.F64[0] = source->peak->xf;
    83         coords->data.F64[1] = source->peak->yf;
     97        float xSource, ySource;         // Coordinates of source
     98        coordsFromSource(&xSource, &ySource, source);
     99
     100        coords->data.F64[0] = xSource;
     101        coords->data.F64[1] = ySource;
    84102
    85103        long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
    … …  
    120138        }
    121139
    122         int x = source->peak->xf + 0.5, y = source->peak->yf + 0.5; // Coordinates of source
     140        float xSource, ySource;         // Coordinates of source
     141        coordsFromSource(&xSource, &ySource, source);
     142
     143        int x = xSource + 0.5, y = ySource + 0.5; // Coordinates of source
    123144
    124145        int xMin = PS_MAX(x - BG_SIZE, 0), xMax = PS_MIN(x + BG_SIZE, numCols);
    … …  
    319340                pmSource *source = pmSourceAlloc();     // Source
    320341                sources->data[0] = source;
    321                 source->peak = pmPeakAlloc(50, 50, 10000, PM_PEAK_LONE);
     342                source->peak = pmPeakAlloc(500, 500, 10000, PM_PEAK_LONE);
    322343                source->psfMag = -13.0;
    323                 pmReadoutFakeFromSources(test, 100, 100, sources, NULL, NULL, psf, 0.1, 0, false, true);
     344                pmReadoutFakeFromSources(test, 1000, 1000, sources, NULL, NULL, psf, 0.1, 0, false, true);
    324345                float sum = 0.0;
    325346                for (int y = 0; y < test->image->numRows; y++) {
    … …  
    351372            }
    352373
     374
     375#if 0
    353376            // Add the background into the target image
    354377            psImage *bgImage = stackBackgroundModel(readout, maskVal, stampSources, footprint); // Image of background
    355378            psBinaryOp(fake->image, fake->image, "+", bgImage);
    356379            psFree(bgImage);
     380#endif
     381
    357382
    358383#ifdef TESTING
    359384            {
     385                pmHDU *hdu = pmHDUFromCell(readout->parent);
    360386                psString name = NULL;
    361387                psStringAppend(&name, "fake_%03d.fits", numInput);
    362388                psFits *fits = psFitsOpen(name, "w");
    363389                psFree(name);
    364                 psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
     390                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
    365391                psFitsClose(fits);
    366392            }
    367393            {
     394                pmHDU *hdu = pmHDUFromCell(readout->parent);
    368395                psString name = NULL;
    369396                psStringAppend(&name, "real_%03d.fits", numInput);
    370397                psFits *fits = psFitsOpen(name, "w");
    371398                psFree(name);
    372                 psFitsWriteImage(fits, NULL, readout->image, 0, NULL);
     399                psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
    373400                psFitsClose(fits);
    374401            }
    … …  
    410437#ifdef TESTING
    411438            {
     439                pmHDU *hdu = pmHDUFromCell(readout->parent);
    412440                psString name = NULL;
    413441                psStringAppend(&name, "conv_%03d.fits", numInput);
    414442                psFits *fits = psFitsOpen(name, "w");
    415443                psFree(name);
    416                 psFitsWriteImage(fits, NULL, output->image, 0, NULL);
     444                psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
    417445                psFitsClose(fits);
    418446            }
    419447            {
     448                pmHDU *hdu = pmHDUFromCell(readout->parent);
    420449                psString name = NULL;
    421450                psStringAppend(&name, "diff_%03d.fits", numInput);
    … …  
    423452                psFree(name);
    424453                psBinaryOp(fake->image, output->image, "-", fake->image);
    425                 psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
     454                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
    426455                psFitsClose(fits);
    427456            }
    … …  
    582611#define RADIUS 10                       // Radius of photometry
    583612#define MIN_ERR 0.05                    // Minimum photometric error, mag
     613#define MAX_MAG -13                     // Maximum magnitude for source
    584614
    585615    // Ensure the normalisation is correct
    … …  
    595625        pmSource *source = sources->data[i]; // Source of interest
    596626        if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) ||
    597             source->errMag > MIN_ERR) {
     627            source->errMag > MIN_ERR || source->psfMag > MAX_MAG) {
    598628            continue;
    599629        }
    600630
    601         float xSrc = source->peak->xf, ySrc = source->peak->yf; // Source coordinates
     631        float xSrc, ySrc;              // Source coordinates
     632        coordsFromSource(&xSrc, &ySrc, source);
    602633        int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x
    603634        int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y
    … …  
    635666#ifdef TESTING
    636667    {
     668        pmHDU *hdu = pmHDUFromCell(readout->parent);
    637669        psString name = NULL;
    638670        psStringAppend(&name, "convolved_%03d.fits", numInput);
    639671        psFits *fits = psFitsOpen(name, "w");
    640672        psFree(name);
    641         psFitsWriteImage(fits, NULL, output->image, 0, NULL);
     673        psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
    642674        psFitsClose(fits);
    643675    }

  • branches/pap_branch_20090108/ppStack/src/ppStackSources.c

    r21094 r21106  
    99//#define TESTING                         // Enable debugging output
    1010
     11// Size of fake image; set by hand because it's trouble to get it from other places
     12#define FAKE_COLS 4861
     13#define FAKE_ROWS 4913
     14
     15
     16
     17#define MIN_MAG -16                     // Minimum magnitude to consider
    1118
    1219float ppStackSourcesTransparency(const psArray *sourceLists, const pmFPAview *view, const pmConfig *config)
    … …  
    5966        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
    6067        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
     68
     69#if 1
     70        // Excise bright sources that might be saturated
     71        psArray *sources = sourceLists->data[i]; // Sources for image
     72        for (int j = 0; j < sources->n; j++) {
     73            pmSource *source = sources->data[j]; // Source of interest
     74            if (source && source->psfMag < MIN_MAG) {
     75                psFree(source);
     76                sources->data[j] = NULL;
     77            }
     78        }
     79#endif
     80
     81
     82#if 1
     83        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
     84        pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
     85        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sourceLists->data[i], NULL, NULL, psf, 5, 0, false, true);
     86        psString name = NULL;
     87        psStringAppend(&name, "start_%03d.fits", i);
     88        psFits *fits = psFitsOpen(name, "w");
     89        psFree(name);
     90        psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
     91        psFitsClose(fits);
     92        psFree(fake);
     93#endif
     94
    6195
    6296        float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time
    … …  
    105139        FILE *outMatches = fopen("source_match.dat", "w"); // Output matches
    106140        psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star
     141        psVector *magErr = psVectorAlloc(num, PS_TYPE_F32); // Errors for each star
    107142        for (int i = 0; i < matches->n; i++) {
    108143            pmSourceMatch *match = matches->data[i]; // Match of interest
    109144            psVectorInit(mag, NAN);
     145            psVectorInit(magErr, NAN);
    110146            for (int j = 0; j < match->num; j++) {
    111147                if (match->mask->data.PS_TYPE_MASK_DATA[j]) {
    … …  
    114150                int index = match->image->data.U32[j]; // Image index
    115151                mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index] - trans->data.F32[index];
     152                magErr->data.F32[index] = match->magErr->data.F32[j];
    116153            }
    117154            for (int j = 0; j < num; j++) {
    118                 fprintf(outMatches, "%f ", mag->data.F32[j]);
     155                fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]);
    119156            }
    120157            fprintf(outMatches, "\n");
    121158        }
    122159        psFree(mag);
     160        psFree(magErr);
    123161        fclose(outMatches);
    124162    }
    … …  
    146184                continue;
    147185            }
     186
     187#if 1
     188            float x = source->modelPSF->params->data.F32[PM_PAR_XPOS],
     189                y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; // Coordinates of source
     190            if (x > 2187 && x < 2191 && y > 2280 && y < 2284) {
     191                fprintf(stderr, "Image %d source: (%f,%f) %f --> %f\n",
     192                        i, x, y, source->psfMag, source->psfMag + magCorr);
     193            }
     194#endif
     195
    148196            source->psfMag += magCorr;
    149197        }
     198
     199#if 1
     200// Size of fake image; set by hand because it's trouble to get it from other places
     201        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
     202        pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
     203        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sources, NULL, NULL, psf, 5, 0, false, true);
     204        psString name = NULL;
     205        psStringAppend(&name, "corrected_%03d.fits", i);
     206        psFits *fits = psFitsOpen(name, "w");
     207        psFree(name);
     208        psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
     209        psFitsClose(fits);
     210        psFree(fake);
     211#endif
     212
    150213    }
    151214    psFree(trans);
    … …  
    162225        psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
    163226                                               transThresh, starRej, starSys);
     227        for (int i = 0; i < num; i++) {
     228            fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
     229        }
    164230        psFree(trans);
    165231    }
Note: See TracChangeset for help on using the changeset viewer.