IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 12, 2008, 4:00:25 PM (18 years ago)
Author:
eugene
Message:

choose the psf 2D scale by measuring the scatter of the E0,1,2 terms, using a bootstrap technique; fill in missing pixels in the psf 2D image

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmPSFtry.c

    r19906 r20078  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2008-10-06 13:05:13 $
     7 *  @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2008-10-13 02:00:22 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3535#include "pmSourcePhotometry.h"
    3636
     37bool printTrendMap (pmTrend2D *trend) {
     38
     39    if (!trend->map) return false;
     40    if (!trend->map->map) return false;
     41
     42    for (int j = 0; j < trend->map->map->numRows; j++) {
     43        for (int i = 0; i < trend->map->map->numCols; i++) {
     44            fprintf (stderr, "%5.2f  ", trend->map->map->data.F32[j][i]);
     45        }
     46        fprintf (stderr, "\t\t\t");
     47        for (int i = 0; i < trend->map->map->numCols; i++) {
     48            fprintf (stderr, "%5.2f  ", trend->map->error->data.F32[j][i]);
     49        }
     50        fprintf (stderr, "\n");
     51    }
     52    return true;
     53}
     54
     55bool psImageMapCleanup (psImageMap *map) {
     56
     57    if ((map->map->numRows == 1) && (map->map->numCols == 1)) return true;
     58
     59    // find the weighted average of all pixels
     60    float Sum = 0.0;
     61    float Wt = 0.0;
     62    for (int j = 0; j < map->map->numRows; j++) {
     63        for (int i = 0; i < map->map->numCols; i++) {
     64            if (!isfinite(map->map->data.F32[j][i])) continue;
     65            Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i];
     66            Wt += map->error->data.F32[j][i];
     67        }
     68    }
     69
     70    float Mean = Sum / Wt;
     71
     72    // do any of the pixels in the map need to be repaired?
     73    // XXX for now, we are just replacing bad pixels with the Mean
     74    for (int j = 0; j < map->map->numRows; j++) {
     75        for (int i = 0; i < map->map->numCols; i++) {
     76            if (isfinite(map->map->data.F32[j][i]) &&
     77                (map->error->data.F32[j][i] > 0.0)) continue;
     78            map->map->data.F32[j][i] = Mean;
     79        }
     80    }
     81    return true;
     82}
     83
    3784// ********  pmPSFtry functions  **************************************************
    3885// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
    … …  
    120167    maskVal |= markVal;
    121168
    122     // stage 1:  fit an EXT model to all candidates PSF sources
     169    // stage 1:  fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF
    123170    psTimerStart ("fit");
    124171    for (int i = 0; i < psfTry->sources->n; i++) {
    … …  
    158205    }
    159206
    160     // stage 2: construct a psf (pmPSF) from this collection of model fits
     207    // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation
    161208    if (!pmPSFFromPSFtry (psfTry)) {
    162209        psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
    … …  
    558605
    559606    if (psf->psfTrendMode == PM_TREND_MAP) {
    560         float errorFloor = 0.0;
    561         float errorTotal = 0.0;
    562         float errorTotalMin = FLT_MAX;
     607        float scatterTotal = 0.0;
     608        float scatterTotalMin = FLT_MAX;
    563609        int entryMin = -1;
    564610
    … …  
    571617
    572618            psVectorInit (mask, 0);
    573             if (!pmPSFFitShapeParamsMap (psf, i, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) {
     619            if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    574620                break;
    575621            }
    576622
    577             // we do not have a good model for the error distributions on e0, e1, e2.
    578             // use the bright end scatter from the constant fit to set the scale for the
    579             // versions with 2D variation.  potentially scale by poisson errors...
    580             // if (i == 1) {
    581             // XXX let's drop this for the moment: relies on valid result from errorFloor
    582             if (0) {
    583                 // if (psf->poissonErrorsParams) do something else..
    584                 dz = psVectorAlloc (sources->n, PS_TYPE_F32);
    585                 for (int i = 0; i < sources->n; i++) {
    586                     dz->data.F32[i] = errorFloor;
    587                 }
    588             }
    589 
    590             // store the resulting errorTotal values and the scales, redo the best
    591             if (errorTotal < errorTotalMin) {
    592                 errorTotalMin = errorTotal;
     623            // store the resulting scatterTotal values and the scales, redo the best
     624            if (scatterTotal < scatterTotalMin) {
     625                scatterTotalMin = scatterTotal;
    593626                entryMin = i;
    594627            }
    … …  
    601634        // XXX supply the resulting mask values back to srcMask
    602635        psVectorInit (mask, 0);
    603         if (!pmPSFFitShapeParamsMap (psf, entryMin, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) {
     636        if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    604637            psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
    605638        }
    … …  
    689722
    690723// fit the shape variations as a psImageMap for the given scale factor
    691 bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *errorFloor, float *errorTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) {
     724bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) {
    692725
    693726    int Nx, Ny;
    694727
    695     // set the map scale to match the aspect ratio
     728    // set the map scale to match the aspect ratio : for a scale of 1, we guarantee
     729    // that we have a single cell
    696730    if (psf->fieldNx > psf->fieldNy) {
    697731        Nx = scale;
    … …  
    707741
    708742    // do we have enough sources for this fine of a grid?
    709     if (x->n < 3*Nx*Ny) {
     743    if (x->n < 10*Nx*Ny) {
    710744        return false;
    711745    }
    712 
    713     // the mask marks the values not used to calculate the ApTrend
    714     psVectorInit (mask, 0);
    715746
    716747    // XXX check this against the exising type
    … …  
    736767    psFree (binning);
    737768
     769    // if the map is 1x1 (a single value), we measure the resulting ensemble scatter
     770
     771    // if the map is more finely sampled, divide the values into two sets: measure the fit from
     772    // one set and the scatter from the other set.
     773    psVector *x_fit = NULL;
     774    psVector *y_fit = NULL;
     775    psVector *x_tst = NULL;
     776    psVector *y_tst = NULL;
     777
     778    psVector *e0obs_fit = NULL;
     779    psVector *e1obs_fit = NULL;
     780    psVector *e2obs_fit = NULL;
     781    psVector *e0obs_tst = NULL;
     782    psVector *e1obs_tst = NULL;
     783    psVector *e2obs_tst = NULL;
     784
     785    if (scale == 1) {
     786        x_fit  = psMemIncrRefCounter (x);
     787        y_fit  = psMemIncrRefCounter (y);
     788        x_tst  = psMemIncrRefCounter (x);
     789        y_tst  = psMemIncrRefCounter (y);
     790        e0obs_fit = psMemIncrRefCounter (e0obs);
     791        e1obs_fit = psMemIncrRefCounter (e1obs);
     792        e2obs_fit = psMemIncrRefCounter (e2obs);
     793        e0obs_tst = psMemIncrRefCounter (e0obs);
     794        e1obs_tst = psMemIncrRefCounter (e1obs);
     795        e2obs_tst = psMemIncrRefCounter (e2obs);
     796    } else {
     797        x_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     798        y_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     799        x_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     800        y_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     801        e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     802        e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     803        e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     804        e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     805        e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     806        e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
     807        for (int i = 0; i < e0obs_fit->n; i++) {
     808            // e0obs->n ==  8 or 9:
     809            // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6
     810            // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7
     811            x_fit->data.F32[i] = x->data.F32[2*i+0]; 
     812            x_tst->data.F32[i] = x->data.F32[2*i+1]; 
     813            y_fit->data.F32[i] = y->data.F32[2*i+0]; 
     814            y_tst->data.F32[i] = y->data.F32[2*i+1]; 
     815
     816            e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0]; 
     817            e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1]; 
     818            e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0];
     819            e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1];
     820            e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0];
     821            e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1];
     822        }
     823    }
     824
     825    // the mask marks the values not used to calculate the ApTrend
     826    psVector *fitMask = psVectorAlloc (x_fit->n, PS_TYPE_U8);
     827    psVectorInit (fitMask, 0);
     828
    738829    // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    739830    // This way, the parameters masked by one of the fits will be applied to the others
    740831    for (int i = 0; i < nIter; i++) {
    741832        // XXX we are using the same stats structure on each pass: do we need to re-init it?
    742         // pmTrend2DFit (psf->params->data[PM_PAR_E0], mask, 0xff, x, y, e0obs, dz);
    743         // psTrace ("psModules.objects", 5, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0obs->n);
    744         // pmTrend2DFit (psf->params->data[PM_PAR_E1], mask, 0xff, x, y, e1obs, dz);
    745         // psTrace ("psModules.objects", 5, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1obs->n);
    746         // pmTrend2DFit (psf->params->data[PM_PAR_E2], mask, 0xff, x, y, e2obs, dz);
    747         // psTrace ("psModules.objects", 5, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2obs->n);
    748 
    749833        psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    750834        psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    … …  
    757841
    758842        trend = psf->params->data[PM_PAR_E0];
    759         status &= pmTrend2DFit (trend, mask, 0xff, x, y, e0obs, dz);
     843        status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL);
    760844        mean = psStatsGetValue (trend->stats, meanOption);
    761845        stdev = psStatsGetValue (trend->stats, stdevOption);
    762         psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs->n);
     846        psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n);
     847        // printTrendMap (trend);
     848        psImageMapCleanup (trend->map);
     849        // printTrendMap (trend);
    763850
    764851        trend = psf->params->data[PM_PAR_E1];
    765         status &= pmTrend2DFit (trend, mask, 0xff, x, y, e1obs, dz);
     852        status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL);
    766853        mean = psStatsGetValue (trend->stats, meanOption);
    767854        stdev = psStatsGetValue (trend->stats, stdevOption);
    768         psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs->n);
     855        psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n);
     856        // printTrendMap (trend);
     857        psImageMapCleanup (trend->map);
     858        // printTrendMap (trend);
    769859
    770860        trend = psf->params->data[PM_PAR_E2];
    771         status &= pmTrend2DFit (trend, mask, 0xff, x, y, e2obs, dz);
     861        status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL);
    772862        mean = psStatsGetValue (trend->stats, meanOption);
    773863        stdev = psStatsGetValue (trend->stats, stdevOption);
    774864        psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
    775 
     865        // printTrendMap (trend);
     866        psImageMapCleanup (trend->map);
     867        // printTrendMap (trend);
    776868    }
    777869    psf->psfTrendStats->clipIter = nIter; // restore default setting
    778870
    779871    // construct the fitted values and the residuals
    780     psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], x, y);
    781     psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], x, y);
    782     psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], x, y);
    783 
    784     psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs, "-", (void *) e0fit);
    785     psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs, "-", (void *) e1fit);
    786     psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs, "-", (void *) e2fit);
    787 
    788 // XXX this code determines the formal error on the map, not the scatter
    789 # if (0)
    790     // measure errors on the map pixels
    791     pmTrend2D *trend;
    792     psStatsOptions meanOption;
    793     psStatsOptions stdevOption;
    794     float mapErrorSum = 0.0;
    795     float mapError = 0.0;
    796 
    797     trend = psf->params->data[PM_PAR_E0];
    798     meanOption = psStatsMeanOption (trend->stats->options);
    799     stdevOption = psStatsStdevOption (trend->stats->options);
    800     mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption));
    801     mapErrorSum += mapError;
    802     psTrace ("psModules.objects", 4, "E0 error: %f\n", sqrt(mapError));
    803 
    804     trend = psf->params->data[PM_PAR_E1];
    805     meanOption = psStatsMeanOption (trend->stats->options);
    806     stdevOption = psStatsStdevOption (trend->stats->options);
    807     mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption));
    808     mapErrorSum += mapError;
    809     psTrace ("psModules.objects", 4, "E1 error: %f\n", sqrt(mapError));
    810 
    811     trend = psf->params->data[PM_PAR_E2];
    812     meanOption = psStatsMeanOption (trend->stats->options);
    813     stdevOption = psStatsStdevOption (trend->stats->options);
    814     mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption));
    815     mapErrorSum += mapError;
    816     psTrace ("psModules.objects", 4, "E2 error: %f\n", sqrt(mapError));
    817 
    818     trend = psf->params->data[PM_PAR_E0];
    819     float mapErrorSum = 0.0;
    820     float mapError = 0.0;
    821     for (int iy = 0; iy < trend->map->error->numRows; iy++) {
    822         for (int ix = 0; ix < trend->map->error->numCols; ix++) {
    823             mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
    824         }
    825     }
    826     psTrace ("psModules.objects", 4, "E0 error: %f\n", sqrt(mapError));
    827     mapErrorSum += mapError;
    828 
    829     trend = psf->params->data[PM_PAR_E1];
    830     mapError = 0.0;
    831     for (int iy = 0; iy < trend->map->error->numRows; iy++) {
    832         for (int ix = 0; ix < trend->map->error->numCols; ix++) {
    833             mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
    834         }
    835     }
    836     psTrace ("psModules.objects", 4, "E1 error: %f\n", sqrt(mapError));
    837     mapErrorSum += mapError;
    838 
    839     trend = psf->params->data[PM_PAR_E2];
    840     mapError = 0.0;
    841     for (int iy = 0; iy < trend->map->error->numRows; iy++) {
    842         for (int ix = 0; ix < trend->map->error->numCols; ix++) {
    843             mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
    844         }
    845     }
    846     psTrace ("psModules.objects", 4, "E2 error: %f\n", sqrt(mapError));
    847     mapErrorSum += mapError;
    848 
    849     psTrace ("psModules.objects", PS_LOG_INFO, "total map error: %f\n", sqrt(mapErrorSum));
    850 # endif
    851 
    852     // measure systematic errorFloor & systematic / photon scale factor
    853     // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
    854     int nGroup = PS_MAX (3*Nx*Ny, 10);
    855     pmPSFShapeParamsErrors (errorFloor, mag, e0res, e1res, e2res, mask, nGroup,
    856                             psStatsStdevOption(psf->psfTrendStats->options));
    857 
    858     // XXX Bogus: *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum);
    859     *errorTotal = *errorFloor;
    860 
    861     psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
    862     psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter floor: %f, error map total: %f\n", *errorFloor, *errorTotal);
     872    psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], x_tst, y_tst);
     873    psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], x_tst, y_tst);
     874    psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], x_tst, y_tst);
     875
     876    psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs_tst, "-", (void *) e0fit);
     877    psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs_tst, "-", (void *) e1fit);
     878    psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs_tst, "-", (void *) e2fit);
     879
     880    // measure scatter for the unfitted points
     881    // psTraceSetLevel ("psLib.math.vectorSampleStdev", 10);
     882    // psTraceSetLevel ("psLib.math.vectorClippedStats", 10);
     883    pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, psStatsStdevOption(psf->psfTrendStats->options));
     884    // psTraceSetLevel ("psLib.math.vectorSampleStdev", 0);
     885    // psTraceSetLevel ("psLib.math.vectorClippedStats", 0);
     886
     887    psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny);
     888    psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter: %f\n", *scatterTotal);
     889
     890    psFree (x_fit);
     891    psFree (y_fit);
     892    psFree (x_tst);
     893    psFree (y_tst);
     894
     895    psFree (e0obs_fit);
     896    psFree (e1obs_fit);
     897    psFree (e2obs_fit);
     898    psFree (e0obs_tst);
     899    psFree (e1obs_tst);
     900    psFree (e2obs_tst);
    863901
    864902    psFree (e0fit);
    … …  
    870908    psFree (e2res);
    871909
     910    psFree (fitMask);
     911
     912    return true;
     913}
     914
     915// calculate the scatter of the parameters
     916bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psStatsOptions stdevOpt)
     917{
     918
     919    // psStats *stats = psStatsAlloc(stdevOpt);
     920    psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_STDEV);
     921
     922    // calculate the root-mean-square of E0, E1, E2
     923    float dEsquare = 0.0;
     924    psStatsInit (stats);
     925    psVectorStats (stats, e0res, NULL, NULL, 0xff);
     926    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
     927
     928    psStatsInit (stats);
     929    psVectorStats (stats, e1res, NULL, NULL, 0xff);
     930    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
     931
     932    psStatsInit (stats);
     933    psVectorStats (stats, e2res, NULL, NULL, 0xff);
     934    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
     935
     936    *scatterTotal = sqrtf(dEsquare);
     937
     938    psFree(stats);
    872939    return true;
    873940}
Note: See TracChangeset for help on using the changeset viewer.