IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 26, 2007, 5:14:57 PM (18 years ago)
Author:
Paul Price
Message:

Adding const in various places.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15562 r15697  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-11-10 01:09:20 $
     7 *  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2007-11-27 03:14:57 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    5656
    5757// allocate a pmPSFtry based on the desired sources and the model (identified by name)
    58 pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options)
     58pmPSFtry *pmPSFtryAlloc (const psArray *sources, const pmPSFOptions *options)
    5959{
    6060    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
    … …  
    9898
    9999// generate a pmPSFtry with a copy of the test PSF sources
    100 pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, pmPSFOptions *options, psMaskType maskVal, psMaskType mark)
     100pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psMaskType maskVal, psMaskType mark)
    101101{
    102102    bool status;
    … …  
    377377
    378378    if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) {
    379         psFree(x);
    380         psFree(y);
    381         psFree(z);
    382         return false;
     379        psFree(x);
     380        psFree(y);
     381        psFree(z);
     382        return false;
    383383    }
    384384
    … …  
    406406        }
    407407
    408         psImageBinning *binning = psImageBinningAlloc();
    409         binning->nXruff = psf->trendNx;
    410         binning->nYruff = psf->trendNy;
    411         binning->nXfine = psf->fieldNx;
    412         binning->nYfine = psf->fieldNy;
    413 
    414         if (psf->psfTrendMode == PM_TREND_MAP) {
    415             psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
    416             psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
    417         }
    418 
    419         // free existing trend, re-alloc
    420         psFree (psf->params->data[i]);
    421         psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats);
    422         psFree (binning);
     408        psImageBinning *binning = psImageBinningAlloc();
     409        binning->nXruff = psf->trendNx;
     410        binning->nYruff = psf->trendNy;
     411        binning->nXfine = psf->fieldNx;
     412        binning->nYfine = psf->fieldNy;
     413
     414        if (psf->psfTrendMode == PM_TREND_MAP) {
     415            psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
     416            psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
     417        }
     418
     419        // free existing trend, re-alloc
     420        psFree (psf->params->data[i]);
     421        psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats);
     422        psFree (binning);
    423423
    424424        // fit the collection of measured parameters to the PSF 2D model
    … …  
    487487        pmSource *source = sources->data[i];
    488488        if (source->modelEXT == NULL) continue;
    489         // XXX I am relying on the fact that none of the masked sources
    490         // have modelEXT set here.  perhaps use the value of psfTry->mask instead?
     489        // XXX I am relying on the fact that none of the masked sources
     490        // have modelEXT set here.  perhaps use the value of psfTry->mask instead?
    491491
    492492        psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
    … …  
    496496        e2->data.F32[i] = pol.e2;
    497497
    498         float flux = source->modelEXT->params->data.F32[PM_PAR_I0];
    499         mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0;
     498        float flux = source->modelEXT->params->data.F32[PM_PAR_I0];
     499        mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0;
    500500    }
    501501
    502502    if (psf->psfTrendMode == PM_TREND_MAP) {
    503         float errorFloor = 0.0;
    504         float errorTotal = 0.0;
    505         float errorTotalMin = FLT_MAX;
    506         int entryMin = -1;
    507 
    508         psVector *dz = NULL;
    509         psVector *mask = psVectorAlloc (sources->n, PS_TYPE_U8);
    510 
    511         // check the fit residuals and increase Nx,Ny until the error is minimized
    512         // pmPSFParamTrend increases the number along the longer of x or y
    513         for (int i = 1; i < PS_MAX (psf->trendNx, psf->trendNy); i++) {
    514 
    515             psVectorInit (mask, 0);
    516             if (!pmPSFFitShapeParamsMap (psf, i, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    517                 break;
    518             }
    519 
    520             // we do not have a good model for the error distributions on e0, e1, e2.
    521             // use the bright end scatter from the constant fit to set the scale for the
    522             // versions with 2D variation.  potentially scale by poisson errors...
    523             if (i == 1) {
    524                 // if (psf->poissonErrorsParams) do something else..
    525                 dz = psVectorAlloc (sources->n, PS_TYPE_F32);
    526                 for (int i = 0; i < sources->n; i++) {
    527                     dz->data.F32[i] = errorFloor;
    528                 }
    529             }
    530 
    531             // store the resulting errorTotal values and the scales, redo the best
    532             if (errorTotal < errorTotalMin) {
    533                 errorTotalMin = errorTotal;
    534                 entryMin = i;
    535             }
    536         }
    537         if (entryMin == -1) {
    538             psError (PS_ERR_UNKNOWN, true, "failed to find image map for shape params");
    539             return false;
    540         }
    541 
    542         // XXX supply the resulting mask values back to srcMask
    543         psVectorInit (mask, 0);
    544         if (!pmPSFFitShapeParamsMap (psf, entryMin, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    545             psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
    546         }
    547        
    548         pmTrend2D *trend = psf->params->data[PM_PAR_E0];
    549         psf->trendNx = trend->map->map->numCols;
    550         psf->trendNy = trend->map->map->numRows;
    551 
    552         psFree (mask);
    553         psFree (dz);
     503        float errorFloor = 0.0;
     504        float errorTotal = 0.0;
     505        float errorTotalMin = FLT_MAX;
     506        int entryMin = -1;
     507
     508        psVector *dz = NULL;
     509        psVector *mask = psVectorAlloc (sources->n, PS_TYPE_U8);
     510
     511        // check the fit residuals and increase Nx,Ny until the error is minimized
     512        // pmPSFParamTrend increases the number along the longer of x or y
     513        for (int i = 1; i < PS_MAX (psf->trendNx, psf->trendNy); i++) {
     514
     515            psVectorInit (mask, 0);
     516            if (!pmPSFFitShapeParamsMap (psf, i, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) {
     517                break;
     518            }
     519
     520            // we do not have a good model for the error distributions on e0, e1, e2.
     521            // use the bright end scatter from the constant fit to set the scale for the
     522            // versions with 2D variation.  potentially scale by poisson errors...
     523            if (i == 1) {
     524                // if (psf->poissonErrorsParams) do something else..
     525                dz = psVectorAlloc (sources->n, PS_TYPE_F32);
     526                for (int i = 0; i < sources->n; i++) {
     527                    dz->data.F32[i] = errorFloor;
     528                }
     529            }
     530
     531            // store the resulting errorTotal values and the scales, redo the best
     532            if (errorTotal < errorTotalMin) {
     533                errorTotalMin = errorTotal;
     534                entryMin = i;
     535            }
     536        }
     537        if (entryMin == -1) {
     538            psError (PS_ERR_UNKNOWN, true, "failed to find image map for shape params");
     539            return false;
     540        }
     541
     542        // XXX supply the resulting mask values back to srcMask
     543        psVectorInit (mask, 0);
     544        if (!pmPSFFitShapeParamsMap (psf, entryMin, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) {
     545            psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
     546        }
     547
     548        pmTrend2D *trend = psf->params->data[PM_PAR_E0];
     549        psf->trendNx = trend->map->map->numCols;
     550        psf->trendNy = trend->map->map->numRows;
     551
     552        psFree (mask);
     553        psFree (dz);
    554554    } else {
    555555
    556         // XXX iterate Nx, Ny based on scatter?
    557         // XXX we force the x & y order to be the same
    558         // modify the order to correspond to the actual number of matched stars:
    559         int order = PS_MAX (psf->trendNx, psf->trendNy);
    560         if ((sources->n < 15) && (order >= 3)) order = 2;
    561         if ((sources->n < 11) && (order >= 2)) order = 1;
    562         if ((sources->n <  8) && (order >= 1)) order = 0;
    563         if ((sources->n <  3)) {
    564             psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
    565             return false;
    566         }
    567         psf->trendNx = order;
    568         psf->trendNy = order;
    569 
    570         // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    571         // This way, the parameters masked by one of the fits will be applied to the others
    572         for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
    573 
    574             // XXX we are using the same stats structure on each pass: do we need to re-init it?
    575             bool status = true;
    576             status &= pmTrend2DFit (psf->params->data[PM_PAR_E0], srcMask, 0xff, x, y, e0, NULL);
    577             psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n);
    578             status &= pmTrend2DFit (psf->params->data[PM_PAR_E1], srcMask, 0xff, x, y, e1, NULL);
    579             psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n);
    580             status &= pmTrend2DFit (psf->params->data[PM_PAR_E2], srcMask, 0xff, x, y, e2, NULL);
    581             psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n);
    582 
    583             if (!status) {
    584                 psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
    585                 return false;
    586             }
    587         }
     556        // XXX iterate Nx, Ny based on scatter?
     557        // XXX we force the x & y order to be the same
     558        // modify the order to correspond to the actual number of matched stars:
     559        int order = PS_MAX (psf->trendNx, psf->trendNy);
     560        if ((sources->n < 15) && (order >= 3)) order = 2;
     561        if ((sources->n < 11) && (order >= 2)) order = 1;
     562        if ((sources->n <  8) && (order >= 1)) order = 0;
     563        if ((sources->n <  3)) {
     564            psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
     565            return false;
     566        }
     567        psf->trendNx = order;
     568        psf->trendNy = order;
     569
     570        // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
     571        // This way, the parameters masked by one of the fits will be applied to the others
     572        for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
     573
     574            // XXX we are using the same stats structure on each pass: do we need to re-init it?
     575            bool status = true;
     576            status &= pmTrend2DFit (psf->params->data[PM_PAR_E0], srcMask, 0xff, x, y, e0, NULL);
     577            psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n);
     578            status &= pmTrend2DFit (psf->params->data[PM_PAR_E1], srcMask, 0xff, x, y, e1, NULL);
     579            psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n);
     580            status &= pmTrend2DFit (psf->params->data[PM_PAR_E2], srcMask, 0xff, x, y, e2, NULL);
     581            psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n);
     582
     583            if (!status) {
     584                psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
     585                return false;
     586            }
     587        }
    588588    }
    589589
    … …  
    591591    if (psTraceGetLevel("psModules.objects") >= 4) {
    592592        FILE *f = fopen ("pol.dat", "w");
    593         fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
     593        fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
    594594        for (int i = 0; i < e0->n; i++) {
    595595            fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
    596596                     x->data.F32[i], y->data.F32[i],
    597                      e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], 
     597                     e0->data.F32[i], e1->data.F32[i], e2->data.F32[i],
    598598                     pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
    599599                     pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]),
    600600                     pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]),
    601                      srcMask->data.U8[i]);
     601                     srcMask->data.U8[i]);
    602602        }
    603603        fclose (f);
    … …  
    615615
    616616    int Nx, Ny;
    617        
     617
    618618    // set the map scale to match the aspect ratio
    619619    if (psf->fieldNx > psf->fieldNy) {
    620         Nx = scale;
    621         Ny = PS_MAX (1, Nx * (psf->fieldNy / psf->fieldNx));
     620        Nx = scale;
     621        Ny = PS_MAX (1, Nx * (psf->fieldNy / psf->fieldNx));
    622622    } else {
    623         Ny = scale;
    624         Nx = PS_MAX (1, Ny * (psf->fieldNx / psf->fieldNy));
     623        Ny = scale;
     624        Nx = PS_MAX (1, Ny * (psf->fieldNx / psf->fieldNy));
    625625    }
    626626
    627627    // do we have enough sources for this fine of a grid?
    628628    if (x->n < 3*Nx*Ny) {
    629         return false;
     629        return false;
    630630    }
    631631
    … …  
    658658    // This way, the parameters masked by one of the fits will be applied to the others
    659659    for (int i = 0; i < nIter; i++) {
    660         // XXX we are using the same stats structure on each pass: do we need to re-init it?
     660        // XXX we are using the same stats structure on each pass: do we need to re-init it?
    661661        pmTrend2DFit (psf->params->data[PM_PAR_E0], mask, 0xff, x, y, e0obs, dz);
    662662        psTrace ("psModules.objects", 5, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0obs->n);
    … …  
    683683    float mapErrorSum = 0.0;
    684684    float mapError = 0.0;
    685     for (int iy = 0; iy < trend->map->error->numRows; iy++) { 
    686         for (int ix = 0; ix < trend->map->error->numCols; ix++) {
    687             mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
    688         }
     685    for (int iy = 0; iy < trend->map->error->numRows; iy++) {
     686        for (int ix = 0; ix < trend->map->error->numCols; ix++) {
     687            mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
     688        }
    689689    }
    690690    psTrace ("psModules.objects", 4, "E0 error: %f\n", sqrt(mapError));
    … …  
    693693    trend = psf->params->data[PM_PAR_E1];
    694694    mapError = 0.0;
    695     for (int iy = 0; iy < trend->map->error->numRows; iy++) { 
    696         for (int ix = 0; ix < trend->map->error->numCols; ix++) {
    697             mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
    698         }
     695    for (int iy = 0; iy < trend->map->error->numRows; iy++) {
     696        for (int ix = 0; ix < trend->map->error->numCols; ix++) {
     697            mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
     698        }
    699699    }
    700700    psTrace ("psModules.objects", 4, "E1 error: %f\n", sqrt(mapError));
    … …  
    703703    trend = psf->params->data[PM_PAR_E2];
    704704    mapError = 0.0;
    705     for (int iy = 0; iy < trend->map->error->numRows; iy++) { 
    706         for (int ix = 0; ix < trend->map->error->numCols; ix++) {
    707             mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
    708         }
     705    for (int iy = 0; iy < trend->map->error->numRows; iy++) {
     706        for (int ix = 0; ix < trend->map->error->numCols; ix++) {
     707            mapError += PS_SQR (trend->map->error->data.F32[iy][ix]);
     708        }
    709709    }
    710710    psTrace ("psModules.objects", 4, "E2 error: %f\n", sqrt(mapError));
    … …  
    742742    int nBin = PS_MAX (mag->n / nGroup, 1);
    743743
    744     // output vectors for ApResid trend 
     744    // output vectors for ApResid trend
    745745    psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32);
    746746
    … …  
    756756    int n = 0;
    757757    for (int i = 0; i < dSo->n; i++) {
    758         int j;
    759         for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
    760             int N = index->data.U32[n];
    761             dE0subset->data.F32[j] = e0res->data.F32[N];
    762             dE1subset->data.F32[j] = e1res->data.F32[N];
    763             dE2subset->data.F32[j] = e2res->data.F32[N];
    764 
    765             mkSubset->data.U8[j]   = mask->data.U8[N];
    766         }
    767         dE0subset->n = j;
    768         dE1subset->n = j;
    769         dE2subset->n = j;
    770         mkSubset->n = j;
    771 
    772         // calculate the root-mean-square of E0, E1, E2
    773         float dEsquare = 0.0;
    774         psStatsInit (statsS);
    775         psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff);
    776         dEsquare += PS_SQR(statsS->robustStdev);
    777 
    778         psStatsInit (statsS);
    779         psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff);
    780         dEsquare += PS_SQR(statsS->robustStdev);
    781 
    782         psStatsInit (statsS);
    783         psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff);
    784         dEsquare += PS_SQR(statsS->robustStdev);
    785 
    786         dSo->data.F32[i] = sqrt(dEsquare);
     758        int j;
     759        for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
     760            int N = index->data.U32[n];
     761            dE0subset->data.F32[j] = e0res->data.F32[N];
     762            dE1subset->data.F32[j] = e1res->data.F32[N];
     763            dE2subset->data.F32[j] = e2res->data.F32[N];
     764
     765            mkSubset->data.U8[j]   = mask->data.U8[N];
     766        }
     767        dE0subset->n = j;
     768        dE1subset->n = j;
     769        dE2subset->n = j;
     770        mkSubset->n = j;
     771
     772        // calculate the root-mean-square of E0, E1, E2
     773        float dEsquare = 0.0;
     774        psStatsInit (statsS);
     775        psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff);
     776        dEsquare += PS_SQR(statsS->robustStdev);
     777
     778        psStatsInit (statsS);
     779        psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff);
     780        dEsquare += PS_SQR(statsS->robustStdev);
     781
     782        psStatsInit (statsS);
     783        psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff);
     784        dEsquare += PS_SQR(statsS->robustStdev);
     785
     786        dSo->data.F32[i] = sqrt(dEsquare);
    787787    }
    788788    psFree (dE0subset);
    … …  
    790790    psFree (dE2subset);
    791791    psFree (mkSubset);
    792    
     792
    793793    psStats *stats = psStatsAlloc (PS_STAT_MIN);
    794794    psVectorStats (stats, dSo, NULL, NULL, 0);
    … …  
    800800
    801801    psFree (dSo);
    802    
     802
    803803    psFree (statsS);
    804804
Note: See TracChangeset for help on using the changeset viewer.