IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 6, 2008, 3:05:13 AM (18 years ago)
Author:
eugene
Message:

upgrades to calculate the PSF clump for subregions on an image; the PSF I/O file now writes the region-based psf clump data, but will read the old format as well

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r19474 r19906  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.61 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2008-09-11 00:38:23 $
     7 *  @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2008-10-06 13:05:13 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    150150    }
    151151    psLogMsg ("psphot.psftry", 4, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("fit"), Next, sources->n);
    152     psTrace ("psphot.psftry", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n);
     152    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n);
    153153
    154154    if (Next == 0) {
    … …  
    215215        psfTry->metricErr->data.F32[i] = source->errMag;
    216216
    217         psTrace ("psphot.psftry", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
     217        psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
    218218        Npsf ++;
    219219    }
    … …  
    221221
    222222    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("fit"), Npsf, sources->n);
    223     psTrace ("psphot.psftry", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n);
     223    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n);
    224224
    225225    if (Npsf == 0) {
    … …  
    568568        // check the fit residuals and increase Nx,Ny until the error is minimized
    569569        // pmPSFParamTrend increases the number along the longer of x or y
    570         for (int i = 1; i < PS_MAX (psf->trendNx, psf->trendNy); i++) {
     570        for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) {
    571571
    572572            psVectorInit (mask, 0);
    … …  
    578578            // use the bright end scatter from the constant fit to set the scale for the
    579579            // versions with 2D variation.  potentially scale by poisson errors...
    580             if (i == 1) {
     580            // if (i == 1) {
     581            // XXX let's drop this for the moment: relies on valid result from errorFloor
     582            if (0) {
    581583                // if (psf->poissonErrorsParams) do something else..
    582584                dz = psVectorAlloc (sources->n, PS_TYPE_F32);
    … …  
    629631        for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
    630632
     633            psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
     634            psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
     635
     636            pmTrend2D *trend = NULL;
     637            float mean, stdev;
     638
    631639            // XXX we are using the same stats structure on each pass: do we need to re-init it?
    632640            bool status = true;
    633             status &= pmTrend2DFit (psf->params->data[PM_PAR_E0], srcMask, 0xff, x, y, e0, NULL);
    634             psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n);
    635             status &= pmTrend2DFit (psf->params->data[PM_PAR_E1], srcMask, 0xff, x, y, e1, NULL);
    636             psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n);
    637             status &= pmTrend2DFit (psf->params->data[PM_PAR_E2], srcMask, 0xff, x, y, e2, NULL);
    638             psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n);
     641
     642            trend = psf->params->data[PM_PAR_E0];
     643            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
     644            mean = psStatsGetValue (trend->stats, meanOption);
     645            stdev = psStatsGetValue (trend->stats, stdevOption);
     646            psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
     647
     648            trend = psf->params->data[PM_PAR_E1];
     649            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
     650            mean = psStatsGetValue (trend->stats, meanOption);
     651            stdev = psStatsGetValue (trend->stats, stdevOption);
     652            psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
     653
     654            trend = psf->params->data[PM_PAR_E2];
     655            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
     656            mean = psStatsGetValue (trend->stats, meanOption);
     657            stdev = psStatsGetValue (trend->stats, stdevOption);
     658            psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
    639659
    640660            if (!status) {
    … …  
    676696    if (psf->fieldNx > psf->fieldNy) {
    677697        Nx = scale;
    678         Ny = PS_MAX (1, Nx * (psf->fieldNy / psf->fieldNx));
     698        float AR = psf->fieldNy / (float) psf->fieldNx;
     699        Ny = (int) (Nx * AR + 0.5);
     700        Ny = PS_MAX (1, Ny);
    679701    } else {
    680702        Ny = scale;
    681         Nx = PS_MAX (1, Ny * (psf->fieldNx / psf->fieldNy));
     703        float AR = psf->fieldNx / (float) psf->fieldNy;
     704        Nx = (int) (Ny * AR + 0.5);
     705        Nx = PS_MAX (1, Nx);
    682706    }
    683707
    … …  
    716740    for (int i = 0; i < nIter; i++) {
    717741        // XXX we are using the same stats structure on each pass: do we need to re-init it?
    718         pmTrend2DFit (psf->params->data[PM_PAR_E0], mask, 0xff, x, y, e0obs, dz);
    719         psTrace ("psModules.objects", 5, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0obs->n);
    720         pmTrend2DFit (psf->params->data[PM_PAR_E1], mask, 0xff, x, y, e1obs, dz);
    721         psTrace ("psModules.objects", 5, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1obs->n);
    722         pmTrend2DFit (psf->params->data[PM_PAR_E2], mask, 0xff, x, y, e2obs, dz);
    723         psTrace ("psModules.objects", 5, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2obs->n);
     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
     749        psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
     750        psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
     751
     752        pmTrend2D *trend = NULL;
     753        float mean, stdev;
     754
     755        // XXX we are using the same stats structure on each pass: do we need to re-init it?
     756        bool status = true;
     757
     758        trend = psf->params->data[PM_PAR_E0];
     759        status &= pmTrend2DFit (trend, mask, 0xff, x, y, e0obs, dz);
     760        mean = psStatsGetValue (trend->stats, meanOption);
     761        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);
     763
     764        trend = psf->params->data[PM_PAR_E1];
     765        status &= pmTrend2DFit (trend, mask, 0xff, x, y, e1obs, dz);
     766        mean = psStatsGetValue (trend->stats, meanOption);
     767        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);
     769
     770        trend = psf->params->data[PM_PAR_E2];
     771        status &= pmTrend2DFit (trend, mask, 0xff, x, y, e2obs, dz);
     772        mean = psStatsGetValue (trend->stats, meanOption);
     773        stdev = psStatsGetValue (trend->stats, stdevOption);
     774        psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
     775
    724776    }
    725777    psf->psfTrendStats->clipIter = nIter; // restore default setting
    … …  
    734786    psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs, "-", (void *) e2fit);
    735787
     788// XXX this code determines the formal error on the map, not the scatter
     789# if (0)
    736790    // measure errors on the map pixels
    737791    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));
    738817
    739818    trend = psf->params->data[PM_PAR_E0];
    … …  
    769848
    770849    psTrace ("psModules.objects", PS_LOG_INFO, "total map error: %f\n", sqrt(mapErrorSum));
     850# endif
    771851
    772852    // measure systematic errorFloor & systematic / photon scale factor
    … …  
    776856                            psStatsStdevOption(psf->psfTrendStats->options));
    777857
    778     *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum);
     858    // XXX Bogus: *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum);
     859    *errorTotal = *errorFloor;
    779860
    780861    psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
    … …  
    815896    for (int i = 0; i < nBin; i++) {
    816897        int j;
     898        int nValid = 0;
    817899        for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
    818900            int N = index->data.U32[n];
    … …  
    822904
    823905            mkSubset->data.U8[j]   = mask->data.U8[N];
    824         }
     906            if (!mask->data.U8[N]) nValid ++;
     907        }
     908        if (nValid < 3) continue;
     909
    825910        dE0subset->n = j;
    826911        dE1subset->n = j;
    … …  
    849934        }
    850935    }
     936    *errorFloor = min;
     937
    851938    psFree (dE0subset);
    852939    psFree (dE1subset);
Note: See TracChangeset for help on using the changeset viewer.