IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 5, 2009, 5:03:33 PM (17 years ago)
Author:
Paul Price
Message:

Merging pap_branch_20090128. Resolved a small number of conflicts. Compiles, but not tested in detail.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotBlendFit.c

    r21359 r21366  
    2121    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    2222    if (!status) {
    23         nThreads = 0;
     23        nThreads = 0;
    2424    }
    2525
    … …  
    4444    psphotInitRadiusPSF (recipe, psf->type);
    4545
    46     // starts the timer, sets up the array of fitSets 
     46    // starts the timer, sets up the array of fitSets
    4747    psphotFitInit (nThreads);
    4848
    … …  
    5656    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
    5757
    58     fprintf (stderr, "starting with %ld sources\n", sources->n);
    59 
    6058    for (int i = 0; i < cellGroups->n; i++) {
    6159
    62         psArray *cells = cellGroups->data[i];
    63 
    64         for (int j = 0; j < cells->n; j++) {
    65 
    66             if (nThreads) {
    67                 // allocate a job -- if threads are not defined, this just runs the job
    68                 psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
    69                 psArray *newSources = psArrayAllocEmpty(16);
    70 
    71                 psArrayAdd(job->args, 1, readout);
    72                 psArrayAdd(job->args, 1, recipe);
    73                 psArrayAdd(job->args, 1, cells->data[j]); // sources
    74                 psArrayAdd(job->args, 1, psf);
    75                 psArrayAdd(job->args, 1, newSources); // return for new sources
    76                 psFree (newSources);
    77 
    78                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
    79                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
    80                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
    81                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
    82 
    83                 if (!psThreadJobAddPending(job)) {
    84                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    85                     psFree (job);
    86                     return NULL;
    87                 }
    88                 psFree(job);
    89             } else {
    90                 int nfit = 0;
    91                 int npsf = 0;
    92                 int next = 0;
    93                 int nfail = 0;
    94                 psArray *newSources = psArrayAllocEmpty(16);
    95 
    96                 if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) {
    97                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    98                     return NULL;
    99                 }
    100                 Nfit += nfit;
    101                 Npsf += npsf;
    102                 Next += next;
    103                 Nfail += nfail;
    104 
    105                 // add these back onto sources
    106                 for (int k = 0; k < newSources->n; k++) {
    107                     psArrayAdd (sources, 16, newSources->data[k]);
    108                 }
    109                 psFree (newSources);
    110             }
    111         }
    112 
    113         if (nThreads) {
    114             // wait for the threads to finish and manage results
    115             if (!psThreadPoolWait (false)) {
    116                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    117                 return NULL;
    118             }
    119 
    120             // we have only supplied one type of job, so we can assume the types here
    121             psThreadJob *job = NULL;
    122             while ((job = psThreadJobGetDone()) != NULL) {
    123                 if (job->args->n < 1) {
    124                     fprintf (stderr, "error with job\n");
    125                 } else {
    126                     psScalar *scalar = NULL;
    127                     scalar = job->args->data[5];
    128                     Nfit += scalar->data.S32;
    129                     scalar = job->args->data[6];
    130                     Npsf += scalar->data.S32;
    131                     scalar = job->args->data[7];
    132                     Next += scalar->data.S32;
    133                     scalar = job->args->data[8];
    134                     Nfail += scalar->data.S32;
    135 
    136                     // add these back onto sources
    137                     psArray *newSources = job->args->data[4];
    138                     for (int j = 0; j < newSources->n; j++) {
    139                         psArrayAdd (sources, 16, newSources->data[j]);
    140                     }
    141                 }
    142                 psFree(job);
    143             }
    144         }
     60        psArray *cells = cellGroups->data[i];
     61
     62        for (int j = 0; j < cells->n; j++) {
     63
     64            if (nThreads) {
     65                // allocate a job -- if threads are not defined, this just runs the job
     66                psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
     67                psArray *newSources = psArrayAllocEmpty(16);
     68
     69                psArrayAdd(job->args, 1, readout);
     70                psArrayAdd(job->args, 1, recipe);
     71                psArrayAdd(job->args, 1, cells->data[j]); // sources
     72                psArrayAdd(job->args, 1, psf);
     73                psArrayAdd(job->args, 1, newSources); // return for new sources
     74                psFree (newSources);
     75
     76                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
     77                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
     78                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
     79                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
     80
     81                if (!psThreadJobAddPending(job)) {
     82                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     83                    psFree (job);
     84                    return NULL;
     85                }
     86                psFree(job);
     87            } else {
     88                int nfit = 0;
     89                int npsf = 0;
     90                int next = 0;
     91                int nfail = 0;
     92                psArray *newSources = psArrayAllocEmpty(16);
     93
     94                if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) {
     95                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     96                    return NULL;
     97                }
     98                Nfit += nfit;
     99                Npsf += npsf;
     100                Next += next;
     101                Nfail += nfail;
     102
     103                // add these back onto sources
     104                for (int k = 0; k < newSources->n; k++) {
     105                    psArrayAdd (sources, 16, newSources->data[k]);
     106                }
     107                psFree (newSources);
     108            }
     109        }
     110
     111        if (nThreads) {
     112            // wait for the threads to finish and manage results
     113            if (!psThreadPoolWait (false)) {
     114                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     115                return NULL;
     116            }
     117
     118            // we have only supplied one type of job, so we can assume the types here
     119            psThreadJob *job = NULL;
     120            while ((job = psThreadJobGetDone()) != NULL) {
     121                if (job->args->n < 1) {
     122                    fprintf (stderr, "error with job\n");
     123                } else {
     124                    psScalar *scalar = NULL;
     125                    scalar = job->args->data[5];
     126                    Nfit += scalar->data.S32;
     127                    scalar = job->args->data[6];
     128                    Npsf += scalar->data.S32;
     129                    scalar = job->args->data[7];
     130                    Next += scalar->data.S32;
     131                    scalar = job->args->data[8];
     132                    Nfail += scalar->data.S32;
     133
     134                    // add these back onto sources
     135                    psArray *newSources = job->args->data[4];
     136                    for (int j = 0; j < newSources->n; j++) {
     137                        psArrayAdd (sources, 16, newSources->data[j]);
     138                    }
     139                }
     140                psFree(job);
     141            }
     142        }
    145143    }
    146144    psFree (cellGroups);
    … …  
    207205        if (source->peak->SN < FIT_SN_LIM) continue;
    208206
    209         // exclude sources outside optional analysis region
     207        // exclude sources outside optional analysis region
    210208        if (source->peak->xf < AnalysisRegion.x0) continue;
    211209        if (source->peak->yf < AnalysisRegion.y0) continue;
    … …  
    216214        if (source->modelPSF == NULL) continue;
    217215
    218         // skip sources which are insignificant flux? 
    219         // XXX this is somewhat ad-hoc
     216        // skip sources which are insignificant flux?
     217        // XXX this is somewhat ad-hoc
    220218        if (source->modelPSF->params->data.F32[1] < 0.1) {
    221219            psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
    … …  
    234232        // try fitting PSFs or extended sources depending on source->mode
    235233        // these functions subtract the resulting fitted source
    236         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    237             if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
    238                 source->type = PM_SOURCE_TYPE_EXTENDED;
    239                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
    240                 Next ++;
    241                 continue;
    242             }
    243         } else {
    244             if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
    245                 source->type = PM_SOURCE_TYPE_STAR;
    246                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
    247                 Npsf ++;
    248                 continue;
    249             }
    250         }
     234        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     235            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
     236                source->type = PM_SOURCE_TYPE_EXTENDED;
     237                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
     238                Next ++;
     239                continue;
     240            }
     241        } else {
     242            if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
     243                source->type = PM_SOURCE_TYPE_STAR;
     244                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
     245                Npsf ++;
     246                continue;
     247            }
     248        }
    251249
    252250        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
    … …  
    317315        if (source->peak->SN < FIT_SN_LIM) continue;
    318316
    319         // exclude sources outside optional analysis region
     317        // exclude sources outside optional analysis region
    320318        if (source->peak->xf < AnalysisRegion.x0) continue;
    321319        if (source->peak->yf < AnalysisRegion.y0) continue;
    … …  
    326324        if (source->modelPSF == NULL) continue;
    327325
    328         // skip sources which are insignificant flux? 
    329         // XXX this is somewhat ad-hoc
     326        // skip sources which are insignificant flux?
     327        // XXX this is somewhat ad-hoc
    330328        if (source->modelPSF->params->data.F32[1] < 0.1) {
    331329            psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
    … …  
    344342        // try fitting PSFs or extended sources depending on source->mode
    345343        // these functions subtract the resulting fitted source
    346         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    347             if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
    348                 source->type = PM_SOURCE_TYPE_EXTENDED;
    349                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
    350                 Next ++;
    351                 continue;
    352             }
    353         } else {
    354             if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
    355                 source->type = PM_SOURCE_TYPE_STAR;
    356                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
    357                 Npsf ++;
    358                 continue;
    359             }
    360         }
     344        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     345            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
     346                source->type = PM_SOURCE_TYPE_EXTENDED;
     347                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
     348                Next ++;
     349                continue;
     350            }
     351        } else {
     352            if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
     353                source->type = PM_SOURCE_TYPE_STAR;
     354                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
     355                Npsf ++;
     356                continue;
     357            }
     358        }
    361359
    362360        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
Note: See TracChangeset for help on using the changeset viewer.