IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15002


Ignore:
Timestamp:
Sep 24, 2007, 6:17:29 PM (19 years ago)
Author:
eugene
Message:

added automatic scale adjustment for PSF parameter model

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15000 r15002  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-09-24 21:27:58 $
     7 *  @version $Revision: 1.47 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2007-09-25 04:17:29 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3535#include "pmSourcePhotometry.h"
    3636
     37int testSaveImage (psMetadata *header, psImage *image, char *filename) {
     38
     39    psFits *fits = psFitsOpen (filename, "w");
     40    psFitsWriteImage (fits, NULL, image, 0, NULL);
     41    psFitsClose (fits);
     42    return (TRUE);
     43}
     44
     45bool pmPSFParamTrend (pmPSF *psf, int scale, float *errorFloor, float *errorTotal, psVector *mask, psVector *x, psVector *y, psVector *e0, psVector *e1, psVector *e2, psVector *dz);
     46bool pmPSFErrorFloor (float *errorFloor, psVector *dMag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup);
     47
    3748// ********  pmPSFtry functions  **************************************************
    3849// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
    … …  
    172183        if (!status) {
    173184            psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL;
    174             psTrace ("psModules.pmPSFtry", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
     185            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
    175186            continue;
    176187        }
    … …  
    179190        if (!pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, mark)) {
    180191            psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
    181             psTrace ("psModules.pmPSFtry", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
     192            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
    182193            continue;
    183194        }
    … …  
    352363    psVector *y  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    353364    psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
     365    psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    354366
    355367    psVector *dz = NULL;
    … …  
    366378        x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
    367379        y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
     380        flux->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_I0];
    368381
    369382        // weight by the error on the source flux
    370383        if (dz) {
    371             dz->data.F32[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0];
     384            // dz->data.F32[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0] / source->modelEXT->params->data.F32[PM_PAR_I0];
     385            dz->data.F32[i] = 1.0;
    372386        }
    373387    }
    … …  
    386400    // threshold.
    387401
    388     // XXX this function needs to check the fit residuals and modify Nx,Ny as needed
    389 
    390402    // convert the measured source shape paramters to polarization terms
    391403    psVector *e0   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    … …  
    403415    }
    404416
    405     // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    406     // This way, the parameters masked by one of the fits will be applied to the others
    407     for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
    408 
    409         // XXX we are using the same stats structure on each pass: do we need to re-init it?
    410         pmTrend2DFit (psf->params->data[PM_PAR_E0], psfTry->mask, 0xff, x, y, e0, dz);
    411         psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n);
    412         pmTrend2DFit (psf->params->data[PM_PAR_E1], psfTry->mask, 0xff, x, y, e1, dz);
    413         psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n);
    414         pmTrend2DFit (psf->params->data[PM_PAR_E2], psfTry->mask, 0xff, x, y, e2, dz);
    415         psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n);
     417    pmTrend2D *trend = psf->params->data[PM_PAR_E0];
     418    if (trend->mode == PM_TREND_MAP) {
     419        float errorFloor = 0.0;
     420        float errorTotal = 0.0;
     421        float errorTotalMin = FLT_MAX;
     422        int entryMin = -1;
     423
     424        // XXX do I need to store the initial value of psfTry->mask and reset?
     425
     426        // check the fit residuals and increase Nx,Ny until the error is minimized
     427        for (int i = 1; i < 10; i++) {
     428
     429            if (!pmPSFParamTrend (psf, i, &errorFloor, &errorTotal, psfTry->mask, x, y, e0, e1, e2, dz)) {
     430                break;
     431            }
     432
     433            // we do not have a good model for the error distributions on e0, e1, e2.
     434            // use the bright end scatter from the constant fit to set the scale for the
     435            // versions with 2D variation
     436            if (i == 1) {
     437                for (int i = 0; i < psfTry->sources->n; i++) {
     438                    if (dz) {
     439                        dz->data.F32[i] = errorFloor;
     440                    }
     441                }
     442            }
     443
     444            // store the resulting errorTotal values and the scales, redo the best
     445            if (errorTotal < errorTotalMin) {
     446                errorTotalMin = errorTotal;
     447                entryMin = i;
     448            }
     449        }
     450        if (entryMin == -1) {
     451            psAbort ("failed on ApResid Trend");
     452        }
     453
     454        // XXX catch error condition
     455        pmPSFParamTrend (psf, entryMin, &errorFloor, &errorTotal, psfTry->mask, x, y, e0, e1, e2, dz);
     456    } else {
     457        // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
     458        // This way, the parameters masked by one of the fits will be applied to the others
     459        for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
     460
     461            // XXX we are using the same stats structure on each pass: do we need to re-init it?
     462            pmTrend2DFit (psf->params->data[PM_PAR_E0], psfTry->mask, 0xff, x, y, e0, dz);
     463            psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n);
     464            pmTrend2DFit (psf->params->data[PM_PAR_E1], psfTry->mask, 0xff, x, y, e1, dz);
     465            psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n);
     466            pmTrend2DFit (psf->params->data[PM_PAR_E2], psfTry->mask, 0xff, x, y, e2, dz);
     467            psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n);
     468        }
    416469    }
    417470
    … …  
    421474        fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
    422475        for (int i = 0; i < e0->n; i++) {
    423             fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
    424                      x->data.F32[i], y->data.F32[i],
     476            fprintf (f, "%f %f  %f  :  %f %f %f  : %f %f %f  : %d\n",
     477                     x->data.F32[i], y->data.F32[i], flux->data.F32[i],
    425478                     e0->data.F32[i], e1->data.F32[i], e2->data.F32[i],
    426479                     pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
    … …  
    437490
    438491    // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
     492    // XXX apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters
    439493    for (int i = 0; i < psf->params->n; i++) {
    440494        switch (i) {
    … …  
    459513        }
    460514
     515        pmTrend2D *trendE0 = psf->params->data[PM_PAR_E0];
     516
     517        // re-create the pmTrend2D using the scales defined for E0
     518        psImageBinning *binning = psImageBinningAlloc();
     519        if (trendE0->mode == PM_TREND_MAP) {
     520            binning->nXruff = trend->map->map->numCols;
     521            binning->nYruff = trend->map->map->numRows;
     522        } else {
     523            binning->nXruff = trend->poly->nX;
     524            binning->nYruff = trend->poly->nY;
     525        }
     526        binning->nXfine = psf->fieldNx;
     527        binning->nYfine = psf->fieldNy;
     528        psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
     529        psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
     530
     531        psFree (psf->params->data[i]);
     532        psf->params->data[i] = pmTrend2DNoImageAlloc (trendE0->mode, binning, psf->psfTrendStats);
     533        psFree (binning);
     534
    461535        // fit the collection of measured parameters to the PSF 2D model
    462536        // the mask is carried from previous steps and updated with this operation
    463537        // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams'
    464         if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, e0, NULL)) {
     538        if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) {
    465539            psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
    466540            psFree(x);
    … …  
    504578    psFree (z);
    505579    psFree (dz);
     580    psFree (flux);
    506581    return true;
    507582}
    508583
     584// XXX this function only works for MAP type -- extend to poly...
     585bool pmPSFParamTrend (pmPSF *psf, int scale, float *errorFloor, float *errorTotal, psVector *mask, psVector *x, psVector *y, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) {
     586
     587    int Nx, Ny;
     588       
     589    // set the map scale to match the aspect ratio
     590    if (psf->fieldNx > psf->fieldNy) {
     591        Nx = scale;
     592        Ny = PS_MAX (1, Nx * (psf->fieldNy / psf->fieldNx));
     593    } else {
     594        Ny = scale;
     595        Nx = PS_MAX (1, Ny * (psf->fieldNx / psf->fieldNy));
     596    }
     597
     598    // do we have enough sources for this fine of a grid?
     599    if (x->n < 3*Nx*Ny) {
     600        return false;
     601    }
     602
     603    // the mask marks the values not used to calculate the ApTrend
     604    psVectorInit (mask, 0);
     605
     606    // XXX check this against the exising type
     607    pmTrend2DMode psfTrendMode = PM_TREND_MAP;
     608
     609    psImageBinning *binning = psImageBinningAlloc();
     610    binning->nXruff = Nx;
     611    binning->nYruff = Ny;
     612    binning->nXfine = psf->fieldNx;
     613    binning->nYfine = psf->fieldNy;
     614    psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
     615    psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
     616
     617    psFree (psf->params->data[PM_PAR_E0]);
     618    psFree (psf->params->data[PM_PAR_E1]);
     619    psFree (psf->params->data[PM_PAR_E2]);
     620
     621    int nIter = psf->psfTrendStats->clipIter;
     622    psf->psfTrendStats->clipIter = 1;
     623    psf->params->data[PM_PAR_E0] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
     624    psf->params->data[PM_PAR_E1] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
     625    psf->params->data[PM_PAR_E2] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
     626    psFree (binning);
     627
     628    // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
     629    // This way, the parameters masked by one of the fits will be applied to the others
     630    for (int i = 0; i < nIter; i++) {
     631        // XXX we are using the same stats structure on each pass: do we need to re-init it?
     632        pmTrend2DFit (psf->params->data[PM_PAR_E0], mask, 0xff, x, y, e0obs, dz);
     633        psTrace ("psModules.objects", 5, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0obs->n);
     634        pmTrend2DFit (psf->params->data[PM_PAR_E1], mask, 0xff, x, y, e1obs, dz);
     635        psTrace ("psModules.objects", 5, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1obs->n);
     636        pmTrend2DFit (psf->params->data[PM_PAR_E2], mask, 0xff, x, y, e2obs, dz);
     637        psTrace ("psModules.objects", 5, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2obs->n);
     638    }
     639    psf->psfTrendStats->clipIter = nIter; // restore default setting
     640
     641    # if (0)
     642    pmTrend2D *trend;
     643    trend = psf->params->data[PM_PAR_E0];
     644    testSaveImage (NULL, trend->map->map, "e0.fits");
     645    testSaveImage (NULL, trend->map->error, "e0d.fits");
     646    trend = psf->params->data[PM_PAR_E1];
     647    testSaveImage (NULL, trend->map->map, "e1.fits");
     648    testSaveImage (NULL, trend->map->error, "e1d.fits");
     649    trend = psf->params->data[PM_PAR_E2];
     650    testSaveImage (NULL, trend->map->map, "e2.fits");
     651    testSaveImage (NULL, trend->map->error, "e2d.fits");
     652    # endif
     653
     654    // construct the fitted values and the residuals
     655    psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], x, y);
     656    psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], x, y);
     657    psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], x, y);
     658
     659    psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs, "-", (void *) e0fit);
     660    psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs, "-", (void *) e1fit);
     661    psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs, "-", (void *) e2fit);
     662
     663    // measure errors on the map pixels
     664    pmTrend2D *trend;
     665
     666    trend = psf->params->data[PM_PAR_E0];
     667    float mapErrorSum = 0.0;
     668    float mapError = 0.0;
     669    for (int iy = 0; iy < trend->map->error->numRows; iy++) {
     670        for (int ix = 0; ix < trend->map->error->numCols; ix++) {
     671            mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
     672        }
     673    }
     674    psTrace ("psModules.objects", 4, "E0 error: %f\n", sqrt(mapError));
     675    mapErrorSum += mapError;
     676
     677    trend = psf->params->data[PM_PAR_E1];
     678    mapError = 0.0;
     679    for (int iy = 0; iy < trend->map->error->numRows; iy++) {
     680        for (int ix = 0; ix < trend->map->error->numCols; ix++) {
     681            mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
     682        }
     683    }
     684    psTrace ("psModules.objects", 4, "E1 error: %f\n", sqrt(mapError));
     685    mapErrorSum += mapError;
     686
     687    trend = psf->params->data[PM_PAR_E2];
     688    mapError = 0.0;
     689    for (int iy = 0; iy < trend->map->error->numRows; iy++) {
     690        for (int ix = 0; ix < trend->map->error->numCols; ix++) {
     691            mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
     692        }
     693    }
     694    psTrace ("psModules.objects", 4, "E2 error: %f\n", sqrt(mapError));
     695    mapErrorSum += mapError;
     696
     697    psTrace ("psModules.objects", PS_LOG_INFO, "total map error: %f\n", sqrt(mapErrorSum));
     698
     699    // measure systematic errorFloor & systematic / photon scale factor
     700    // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
     701    int nGroup = PS_MAX (3*Nx*Ny, 10);
     702    pmPSFErrorFloor (errorFloor, dz, e0res, e1res, e2res, mask, nGroup);
     703
     704    *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum);
     705
     706    psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
     707    psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter floor: %f, error map total: %f\n", *errorFloor, *errorTotal);
     708
     709    psFree (e0fit);
     710    psFree (e1fit);
     711    psFree (e2fit);
     712
     713    psFree (e0res);
     714    psFree (e1res);
     715    psFree (e2res);
     716
     717    return true;
     718}
     719
     720bool pmPSFErrorFloor (float *errorFloor, psVector *dMag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup) {
     721
     722    psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     723
     724    // measure the trend in bins with 10 values each; if < 10 total, use them all
     725    int nBin = PS_MAX (dMag->n / nGroup, 1);
     726
     727    // output vectors for ApResid trend
     728    psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32);
     729
     730    // use dMag to group the dMag and dap vectors
     731    psVector *index = psVectorSortIndex (NULL, dMag);
     732
     733    // subset vectors for dMag and dap values within the given range
     734    psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
     735    psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
     736    psVector *dE2subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
     737    psVector *mkSubset  = psVectorAllocEmpty (nGroup, PS_TYPE_U8);
     738
     739    int n = 0;
     740    for (int i = 0; i < dSo->n; i++) {
     741        int j;
     742        for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) {
     743            int N = index->data.U32[n];
     744            dE0subset->data.F32[j] = e0res->data.F32[N];
     745            dE1subset->data.F32[j] = e1res->data.F32[N];
     746            dE2subset->data.F32[j] = e2res->data.F32[N];
     747
     748            mkSubset->data.U8[j]   = mask->data.U8[N];
     749        }
     750        dE0subset->n = j;
     751        dE1subset->n = j;
     752        dE2subset->n = j;
     753        mkSubset->n = j;
     754
     755        // calculate the root-mean-square of E0, E1, E2
     756        float dEsquare = 0.0;
     757        psStatsInit (statsS);
     758        psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff);
     759        dEsquare += PS_SQR(statsS->robustStdev);
     760
     761        psStatsInit (statsS);
     762        psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff);
     763        dEsquare += PS_SQR(statsS->robustStdev);
     764
     765        psStatsInit (statsS);
     766        psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff);
     767        dEsquare += PS_SQR(statsS->robustStdev);
     768
     769        dSo->data.F32[i] = sqrt(dEsquare);
     770    }
     771    psFree (dE0subset);
     772    psFree (dE1subset);
     773    psFree (dE2subset);
     774    psFree (mkSubset);
     775   
     776    psStats *stats = psStatsAlloc (PS_STAT_MIN);
     777    psVectorStats (stats, dSo, NULL, NULL, 0);
     778
     779    *errorFloor = stats->min;
     780
     781    psFree (stats);
     782    psFree (index);
     783
     784    psFree (dSo);
     785   
     786    psFree (statsS);
     787
     788    return true;
     789}
Note: See TracChangeset for help on using the changeset viewer.