IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28405 for trunk/psphot


Ignore:
Timestamp:
Jun 18, 2010, 2:25:38 PM (16 years ago)
Author:
Paul Price
Message:

Discovered a race condition in the threading process. The thread process was: psThreadJobAlloc, populate job with arguments, psThreadJobAddPending, psFree the job. This is not thread-safe, because we're decrementing the job's reference counter while the reference count is also being fiddled within the thread system (psThread.c; by pulling the job off the 'pending' queue and putting it on the 'done' queue). The reference count fiddling within the thread system was protected by a mutex, but the external free was not. Therefore, I'm pulling the free into the thread system so it can be protected (without the user having to worry about locks to do the simple case). This means that the user should no longer touch the job once he hands it over to the thread system (unless it's protected by calls to psThreadLock/psThreadUnlock, or it's on the other side of psThreadPoolWait from all the adds so that the threads aren't doing anything).

Location:
trunk/psphot/src
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/psphot/src/psphotApResid.c

    r28013 r28405  
    2222    // loop over the available readouts
    2323    for (int i = 0; i < num; i++) {
    24         if (i == chisqNum) continue; // skip chisq image
    25         if (!psphotApResidReadout (config, view, filerule, i, recipe)) {
     24        if (i == chisqNum) continue; // skip chisq image
     25        if (!psphotApResidReadout (config, view, filerule, i, recipe)) {
    2626            psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for %s entry %d", filerule, i);
    27             return false;
    28         }
     27            return false;
     28        }
    2929    }
    3030    return true;
     
    5656
    5757    if (!sources->n) {
    58         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid");
    59         return true;
     58        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid");
     59        return true;
    6060    }
    6161
     
    6666    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    6767    if (!status) {
    68         nThreads = 0;
     68        nThreads = 0;
    6969    }
    7070
     
    128128    for (int i = 0; i < cellGroups->n; i++) {
    129129
    130         psArray *cells = cellGroups->data[i];
    131 
    132         for (int j = 0; j < cells->n; j++) {
    133 
    134             // allocate a job -- if threads are not defined, this just runs the job
    135             psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS");
    136 
    137             psArrayAdd(job->args, 1, cells->data[j]); // sources
    138             psArrayAdd(job->args, 1, psf);
    139             PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
    140             PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
    141             PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
    142 
    143             PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
    144             PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nfail
    145 
    146             if (!psThreadJobAddPending(job)) {
    147                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    148                 psFree (job);
    149                 return false;
    150             }
    151             psFree(job);
    152         }
    153 
    154         // wait for the threads to finish and manage results
    155         if (!psThreadPoolWait (false)) {
    156             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    157             return false;
    158         }
    159 
    160         // we have only supplied one type of job, so we can assume the types here
    161         psThreadJob *job = NULL;
    162         while ((job = psThreadJobGetDone()) != NULL) {
    163             if (job->args->n < 1) {
    164                 fprintf (stderr, "error with job\n");
    165             } else {
    166                 psScalar *scalar = NULL;
    167                 scalar = job->args->data[5];
    168                 Nskip += scalar->data.S32;
    169                 scalar = job->args->data[6];
    170                 Nfail += scalar->data.S32;
    171             }
    172             psFree(job);
    173         }
     130        psArray *cells = cellGroups->data[i];
     131
     132        for (int j = 0; j < cells->n; j++) {
     133
     134            // allocate a job -- if threads are not defined, this just runs the job
     135            psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS");
     136
     137            psArrayAdd(job->args, 1, cells->data[j]); // sources
     138            psArrayAdd(job->args, 1, psf);
     139            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
     140            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     141            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
     142
     143            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
     144            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nfail
     145
     146            if (!psThreadJobAddPending(job)) {
     147                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     148                return false;
     149            }
     150        }
     151
     152        // wait for the threads to finish and manage results
     153        if (!psThreadPoolWait (false)) {
     154            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     155            return false;
     156        }
     157
     158        // we have only supplied one type of job, so we can assume the types here
     159        psThreadJob *job = NULL;
     160        while ((job = psThreadJobGetDone()) != NULL) {
     161            if (job->args->n < 1) {
     162                fprintf (stderr, "error with job\n");
     163            } else {
     164                psScalar *scalar = NULL;
     165                scalar = job->args->data[5];
     166                Nskip += scalar->data.S32;
     167                scalar = job->args->data[6];
     168                Nfail += scalar->data.S32;
     169            }
     170            psFree(job);
     171        }
    174172    }
    175173
     
    184182    Npsf = 0;
    185183
    186 # ifdef DEBUG   
     184# ifdef DEBUG
    187185    FILE *f = fopen ("apresid.dat", "w");
    188186    psAssert (f, "failed open");
     
    199197        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
    200198
    201         if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
    202         if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
    203         if (source->mode &  PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT");
    204            
     199        if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
     200        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
     201        if (source->mode &  PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT");
     202
    205203        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    206204            continue;
    207205        }
    208206
    209         // XXX make this user-configurable?
    210         if (source->errMag > 0.01) continue;
     207        // XXX make this user-configurable?
     208        if (source->errMag > 0.01) continue;
    211209
    212210        // aperture residual for this source
     
    221219
    222220# ifdef DEBUG
    223         fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
    224                  source->peak->xf, source->peak->yf,
    225                  source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS],
    226                 source->psfMag, source->apMag, source->errMag,
    227                  source->modelPSF->params->data.F32[PM_PAR_I0],
    228                  source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY],
    229                 source->modelPSF->params->data.F32[PM_PAR_7]);
     221        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
     222                 source->peak->xf, source->peak->yf,
     223                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS],
     224                source->psfMag, source->apMag, source->errMag,
     225                 source->modelPSF->params->data.F32[PM_PAR_I0],
     226                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY],
     227                source->modelPSF->params->data.F32[PM_PAR_7]);
    230228# endif
    231         if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
    232         if (!isfinite(source->errMag)) psAbort ("nan in errMag");
    233         if (!isfinite(source->apMag)) psAbort ("nan in apMag");
    234         if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
    235         if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos");
     229        if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
     230        if (!isfinite(source->errMag)) psAbort ("nan in errMag");
     231        if (!isfinite(source->apMag)) psAbort ("nan in apMag");
     232        if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
     233        if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos");
    236234
    237235        psVectorAppend (mag, source->psfMag);
     
    253251    if (Npsf < APTREND_NSTAR_MIN) {
    254252        psWarning("Only %d valid aperture residual sources (need %d), giving up", Npsf, APTREND_NSTAR_MIN);
    255         goto escape;
    256     }
    257 
    258     // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different 
    259     // definition for 'order' (order_MAP = order_POLY + 1).  in addition, we have a 
     253        goto escape;
     254    }
     255
     256    // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different
     257    // definition for 'order' (order_MAP = order_POLY + 1).  in addition, we have a
    260258    // user-specified MAX order, which we should respect, regardless of the mode
    261259
     
    270268    pmTrend2DMode mode = PM_TREND_MAP;
    271269    if (mode == PM_TREND_MAP) {
    272         MaxOrderForStars ++;
    273     } 
     270        MaxOrderForStars ++;
     271    }
    274272    APTREND_ORDER_MAX = PS_MIN (APTREND_ORDER_MAX, MaxOrderForStars);
    275273
     
    283281    int NY = readout->image->numRows;
    284282    for (int i = 1; i <= APTREND_ORDER_MAX; i++) {
    285        
    286         int Nx, Ny;
    287         if (NX > NY) {
    288             Nx = i;
    289             Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5));
    290         } else {
    291             Ny = i;
    292             Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5));
    293         }
    294 
    295         float errorFloor;
     283
     284        int Nx, Ny;
     285        if (NX > NY) {
     286            Nx = i;
     287            Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5));
     288        } else {
     289            Ny = i;
     290            Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5));
     291        }
     292
     293        float errorFloor;
    296294        pmTrend2D *apTrend = psphotApResidTrend (&errorFloor, readout, Nx, Ny, xPos, yPos, apResid, dMag);
    297         if (!apTrend) {
    298             continue;
    299         }
    300 
    301         // apply ApTrend results
    302         // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
    303         // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
    304         // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
    305         // if (!isfinite(ApResid)) psAbort("nan apresid @ center");
     295        if (!apTrend) {
     296            continue;
     297        }
     298
     299        // apply ApTrend results
     300        // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
     301        // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
     302        // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
     303        // if (!isfinite(ApResid)) psAbort("nan apresid @ center");
    306304
    307305        // store the minimum errorFloor and best ApTrend to keep
    308306        if (errorFloor < errorFloorMin) {
    309307            errorFloorMin = errorFloor;
    310             psFree (psf->ApTrend);
    311             psf->ApTrend = psMemIncrRefCounter(apTrend);
    312         }
    313         psFree (apTrend);
     308            psFree (psf->ApTrend);
     309            psf->ApTrend = psMemIncrRefCounter(apTrend);
     310        }
     311        psFree (apTrend);
    314312    }
    315313    if (psf->ApTrend == NULL) {
    316314        psWarning("Failed to find a valid aperture residual value");
    317         goto escape;
     315        goto escape;
    318316    }
    319317
     
    383381    if (!pmTrend2DFit (apTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft)) {
    384382        psWarning("Failed to fit trend for %d x %d map", Nx, Ny);
    385         psFree (apTrend);
    386         return NULL;
     383        psFree (apTrend);
     384        return NULL;
    387385    }
    388386    if (apTrend->mode == PM_TREND_MAP) {
    389         // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST:
    390         psImageMapRepair (apTrend->map->map);
    391         // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST:
     387        // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST:
     388        psImageMapRepair (apTrend->map->map);
     389        // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST:
    392390    }
    393391
     
    400398    if (!isfinite(*apResidSysErr)) {
    401399        psWarning("Failed to find systematic error for %d x %d map", Nx, Ny);
    402         psFree (apTrend);
    403         return NULL;
     400        psFree (apTrend);
     401        return NULL;
    404402    }
    405403
     
    408406
    409407    if (psTraceGetLevel("psphot") >= 4) {
    410         char filename[64];
    411         snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny);
     408        char filename[64];
     409        snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny);
    412410        FILE *dumpFile = fopen (filename, "w");
    413411        for (int i = 0; i < xPos->n; i++) {
     
    457455        }
    458456
    459         // clear the mask bit and set the circular mask pixels
    460         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    461         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
    462 
    463         bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
    464 
    465         // clear the mask bit
    466         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     457        // clear the mask bit and set the circular mask pixels
     458        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     459        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
     460
     461        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
     462
     463        // clear the mask bit
     464        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    467465
    468466        // re-subtract the object, leave local sky
    469467        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    470468
    471         if (!status) {
    472             Nskip ++;
    473             psTrace ("psphot", 3, "skip : bad source mag");
    474             continue;
    475         }
    476    
    477         if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    478             Nfail ++;
    479             psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
    480             continue;
    481         }
    482         source->mode |= PM_SOURCE_MODE_AP_MAGS;
     469        if (!status) {
     470            Nskip ++;
     471            psTrace ("psphot", 3, "skip : bad source mag");
     472            continue;
     473        }
     474
     475        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
     476            Nfail ++;
     477            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
     478            continue;
     479        }
     480        source->mode |= PM_SOURCE_MODE_AP_MAGS;
    483481    }
    484482
  • trunk/psphot/src/psphotBlendFit.c

    r28013 r28405  
    1515    // loop over the available readouts
    1616    for (int i = 0; i < num; i++) {
    17         if (!psphotBlendFitReadout (config, view, filerule, i, recipe)) {
     17        if (!psphotBlendFitReadout (config, view, filerule, i, recipe)) {
    1818            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (non-linear) for %s entry %d", filerule, i);
    19             return false;
    20         }
     19            return false;
     20        }
    2121    }
    2222    return true;
     
    4848
    4949    if (!sources->n) {
    50         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit");
    51         return true;
     50        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit");
     51        return true;
    5252    }
    5353
     
    8787    sources = psArraySort (sources, pmSourceSortBySN);
    8888    if (!sources->n) {
    89         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
    90         return true;
     89        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
     90        return true;
    9191    }
    9292
     
    103103        for (int j = 0; j < cells->n; j++) {
    104104
    105             // allocate a job -- if threads are not defined, this just runs the job
    106             psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
    107             psArray *newSources = psArrayAllocEmpty(16);
    108 
    109             psArrayAdd(job->args, 1, readout);
    110             psArrayAdd(job->args, 1, recipe);
    111             psArrayAdd(job->args, 1, cells->data[j]); // sources
    112             psArrayAdd(job->args, 1, psf);
    113             psArrayAdd(job->args, 1, newSources); // return for new sources
    114             psFree (newSources);
    115 
    116             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
    117             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
    118             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
    119             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
    120 
    121             if (!psThreadJobAddPending(job)) {
    122                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    123                 psFree (job);
    124                 return NULL;
    125             }
    126             psFree(job);
     105            // allocate a job -- if threads are not defined, this just runs the job
     106            psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
     107            psArray *newSources = psArrayAllocEmpty(16);
     108
     109            psArrayAdd(job->args, 1, readout);
     110            psArrayAdd(job->args, 1, recipe);
     111            psArrayAdd(job->args, 1, cells->data[j]); // sources
     112            psArrayAdd(job->args, 1, psf);
     113            psArrayAdd(job->args, 1, newSources); // return for new sources
     114            psFree (newSources);
     115
     116            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
     117            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
     118            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
     119            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
     120
     121            if (!psThreadJobAddPending(job)) {
     122                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     123                return NULL;
     124            }
    127125
    128126# if (0)
     
    152150        }
    153151
    154         // wait for the threads to finish and manage results
    155         if (!psThreadPoolWait (false)) {
    156             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    157             return NULL;
    158         }
    159 
    160         // we have only supplied one type of job, so we can assume the types here
    161         psThreadJob *job = NULL;
    162         while ((job = psThreadJobGetDone()) != NULL) {
    163             if (job->args->n < 1) {
    164                 fprintf (stderr, "error with job\n");
    165             } else {
    166                 psScalar *scalar = NULL;
    167                 scalar = job->args->data[5];
    168                 Nfit += scalar->data.S32;
    169                 scalar = job->args->data[6];
    170                 Npsf += scalar->data.S32;
    171                 scalar = job->args->data[7];
    172                 Next += scalar->data.S32;
    173                 scalar = job->args->data[8];
    174                 Nfail += scalar->data.S32;
    175 
    176                 // add these back onto sources
    177                 psArray *newSources = job->args->data[4];
    178                 for (int j = 0; j < newSources->n; j++) {
    179                     psArrayAdd (sources, 16, newSources->data[j]);
    180                 }
    181             }
    182             psFree(job);
     152        // wait for the threads to finish and manage results
     153        if (!psThreadPoolWait (false)) {
     154            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     155            return NULL;
     156        }
     157
     158        // we have only supplied one type of job, so we can assume the types here
     159        psThreadJob *job = NULL;
     160        while ((job = psThreadJobGetDone()) != NULL) {
     161            if (job->args->n < 1) {
     162                fprintf (stderr, "error with job\n");
     163            } else {
     164                psScalar *scalar = NULL;
     165                scalar = job->args->data[5];
     166                Nfit += scalar->data.S32;
     167                scalar = job->args->data[6];
     168                Npsf += scalar->data.S32;
     169                scalar = job->args->data[7];
     170                Next += scalar->data.S32;
     171                scalar = job->args->data[8];
     172                Nfail += scalar->data.S32;
     173
     174                // add these back onto sources
     175                psArray *newSources = job->args->data[4];
     176                for (int j = 0; j < newSources->n; j++) {
     177                    psArrayAdd (sources, 16, newSources->data[j]);
     178                }
     179            }
     180            psFree(job);
    183181            }
    184182    }
     
    278276                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
    279277                Next ++;
    280                 source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
     278                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    281279                continue;
    282280            }
     
    286284                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
    287285                Npsf ++;
    288                 source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
     286                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    289287                continue;
    290288            }
  • trunk/psphot/src/psphotGuessModels.c

    r28013 r28405  
    2121    // loop over the available readouts
    2222    for (int i = 0; i < num; i++) {
    23         if (i == chisqNum) continue; // skip chisq image
     23        if (i == chisqNum) continue; // skip chisq image
    2424        if (!psphotGuessModelsReadout (config, view, filerule, i)) {
    2525            psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for %s entry %d", filerule, i);
     
    106106            if (!psThreadJobAddPending(job)) {
    107107                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    108                 psFree (job);
    109108                return false;
    110109            }
    111             psFree(job);
    112110        }
    113111
  • trunk/psphot/src/psphotMagnitudes.c

    r28098 r28405  
    1818    // loop over the available readouts
    1919    for (int i = 0; i < num; i++) {
    20         if (i == chisqNum) continue; // skip chisq image
    21 
    22         // find the currently selected readout
    23         pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
    24         psAssert (file, "missing file?");
    25 
    26         pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    27         psAssert (readout, "missing readout?");
    28 
    29         pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    30         psAssert (detections, "missing detections?");
    31 
    32         psArray *sources = detections->allSources;
    33         psAssert (sources, "missing sources?");
    34 
    35         pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    36         psAssert (psf, "missing psf?");
    37 
    38         if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) {
     20        if (i == chisqNum) continue; // skip chisq image
     21
     22        // find the currently selected readout
     23        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
     24        psAssert (file, "missing file?");
     25
     26        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     27        psAssert (readout, "missing readout?");
     28
     29        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     30        psAssert (detections, "missing detections?");
     31
     32        psArray *sources = detections->allSources;
     33        psAssert (sources, "missing sources?");
     34
     35        pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     36        psAssert (psf, "missing psf?");
     37
     38        if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) {
    3939            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
    40             return false;
    41         }
     40            return false;
     41        }
    4242    }
    4343    return true;
     
    5050
    5151    if (!sources->n) {
    52         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source magnitudes");
    53         return true;
    54     }
    55        
     52        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source magnitudes");
     53        return true;
     54    }
     55
    5656    psTimerStart ("psphot.mags");
    5757
     
    122122            if (!psThreadJobAddPending(job)) {
    123123                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    124                 psFree (job);
    125124                return false;
    126125            }
    127             psFree(job);
    128126
    129127# if (0)
     
    184182        }
    185183
    186         // clear the mask bit and set the circular mask pixels
    187         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    188         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
     184        // clear the mask bit and set the circular mask pixels
     185        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     186        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
    189187
    190188        status = pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
    191189        if (status && isfinite(source->apMag)) Nap ++;
    192190
    193         // clear the mask bit
    194         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     191        // clear the mask bit
     192        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    195193
    196194        // re-subtract the object, leave local sky
     
    273271            if (!psThreadJobAddPending(job)) {
    274272                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    275                 psFree (job);
    276273                return false;
    277274            }
    278             psFree(job);
    279 
    280275        }
    281276
  • trunk/psphot/src/psphotSourceStats.c

    r28013 r28405  
    33// convert detections to sources and measure their basic properties (moments, local sky, sky
    44// variance) Note: this function only generates sources for the new peaks (peak->assigned).
    5 // The new sources are added to any existing sources on detections->newSources.  The sources 
     5// The new sources are added to any existing sources on detections->newSources.  The sources
    66// on detections->allSources are ignored.
    77bool psphotSourceStats (pmConfig *config, const pmFPAview *view, const char *filerule, bool setWindow)
     
    2222    // loop over the available readouts
    2323    for (int i = 0; i < num; i++) {
    24         // if (i == chisqNum) continue; // skip chisq image
     24        // if (i == chisqNum) continue; // skip chisq image
    2525        if (!psphotSourceStatsReadout (config, view, filerule, i, recipe, setWindow)) {
    2626            psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for %s entry %d", filerule, i);
     
    9191    // generate the array of sources, define the associated pixel
    9292    if (!detections->newSources) {
    93         detections->newSources = psArrayAllocEmpty (peaks->n);
     93        detections->newSources = psArrayAllocEmpty (peaks->n);
    9494    }
    9595    sources = detections->newSources;
     
    107107        // create a new source
    108108        pmSource *source = pmSourceAlloc();
    109         source->imageID = index;
     109        source->imageID = index;
    110110
    111111        // add the peak
     
    184184            if (!psThreadJobAddPending(job)) {
    185185                psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
    186                 psFree (job);
    187186                psFree(detections->newSources);
    188187                return false;
    189188            }
    190             psFree(job);
    191189        }
    192190
     
    194192        if (!psThreadPoolWait (false)) {
    195193            psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
    196             psFree(detections->newSources);
     194            psFree(detections->newSources);
    197195            return false;
    198196        }
     
    315313            if (!psThreadJobAddPending(job)) {
    316314                psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
    317                 psFree (job);
    318315                return NULL;
    319316            }
    320             psFree(job);
    321317        }
    322318
     
    380376        pmSource *source = sources->data[i];
    381377
    382         if (source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue;
    383         source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
     378        if (source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue;
     379        source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
    384380
    385381        // skip faint sources for moments measurement
  • trunk/psphot/src/psphotTest.c

    r24187 r28405  
    2525
    2626bool FillImage_Threaded (psThreadJob *job) {
    27    
     27
    2828    psImage *image = job->args->data[0];
    2929    int xs = PS_SCALAR_VALUE(job->args->data[1],S32);
     
    3636    psRegion region = psRegionSet (xs, xs + dx, ys, ys + dy);
    3737    for (int i = 0; i < 100; i++) {
    38         psImage *subset = psImageSubset (image, region);
    39         psImageInit (subset, value + i);
    40         psFree (subset);
     38        psImage *subset = psImageSubset (image, region);
     39        psImageInit (subset, value + i);
     40        psFree (subset);
    4141    }
    4242    return true;
     
    4646
    4747    if (argc != 3) {
    48         fprintf (stderr, "USAGE: psphotTest (output.fits) (nThreads)\n");
    49         exit (2);
     48        fprintf (stderr, "USAGE: psphotTest (output.fits) (nThreads)\n");
     49        exit (2);
    5050    }
    5151
     
    6262
    6363    for (int ix = 0; ix < 1000; ix += 100) {
    64         for (int iy = 0; iy < 1000; iy += 100) {
     64        for (int iy = 0; iy < 1000; iy += 100) {
    6565
    66             // allocate a job -- if threads are not defined, this just runs the job
    67             psThreadJob *job = psThreadJobAlloc ("FILL_IMAGE");
     66            // allocate a job -- if threads are not defined, this just runs the job
     67            psThreadJob *job = psThreadJobAlloc ("FILL_IMAGE");
    6868
    69             psArrayAdd(job->args, 1, image);
    70             PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32);
    71             PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32);
    72             PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32);
    73             PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32);
    74             PS_ARRAY_ADD_SCALAR(job->args, ix + iy, PS_TYPE_S32);
     69            psArrayAdd(job->args, 1, image);
     70            PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32);
     71            PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32);
     72            PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32);
     73            PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32);
     74            PS_ARRAY_ADD_SCALAR(job->args, ix + iy, PS_TYPE_S32);
    7575
    76             // FillImage (image, ix, iy, 100, 100, ix + iy);
     76            // FillImage (image, ix, iy, 100, 100, ix + iy);
    7777
    78             if (!psThreadJobAddPending(job)) {
    79                 fprintf (stderr, "failure to run FillImage(1)");
    80                 psFree (job);
    81                 exit (1);
    82             }
    83             psFree(job);
    84         }
     78            if (!psThreadJobAddPending(job)) {
     79                fprintf (stderr, "failure to run FillImage(1)");
     80                exit (1);
     81            }
     82        }
    8583    }
    8684
     
    8886    // wait for the threads to finish and manage results
    8987    if (!psThreadPoolWait (true)) {
    90         fprintf (stderr, "failure to run FillImage (2)");
    91         exit (1);
     88        fprintf (stderr, "failure to run FillImage (2)");
     89        exit (1);
    9290    }
    9391
Note: See TracChangeset for help on using the changeset viewer.