IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 32023


Ignore:
Timestamp:
Aug 3, 2011, 6:29:49 PM (15 years ago)
Author:
eugene
Message:

thread psphotExtendedAnalsis and psphotRadialApertures; subtract the fake source image so output residuals image is clean

Location:
branches/eam_branches/ipp-20110710/psphot/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110710/psphot/src/psphot.h

    r31990 r32023  
    117117bool            psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule);
    118118bool            psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     119bool            psphotExtendedSourceAnalysis_Threaded (psThreadJob *job);
    119120
    120121bool            psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule);
    … …  
    445446bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule);
    446447bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     448bool psphotRadialApertures_Threaded (psThreadJob *job);
    447449bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry);
    448450

  • branches/eam_branches/ipp-20110710/psphot/src/psphotEfficiency.c

    r31452 r32023  
    9999
    100100/// Generate a fake image and add it in to the existing readout
    101 static bool effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources
     101static pmReadout *effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources
    102102                        const pmReadout *ro,            // Readout of interest
    103103                        const pmPSF *psf,               // Point-spread function
    … …  
    152152        psFree(xAll);
    153153        psFree(yAll);
    154         return false;
     154        return NULL;
    155155    }
    156156    psFree(magAll);
    … …  
    161161
    162162    psBinaryOp(ro->image, ro->image, "+", fakeRO->image);
    163     psFree(fakeRO);
    164 
    165     return true;
     163
     164    // return the readout so we can subtract it later
     165    return fakeRO;
    166166}
    167167
    … …  
    290290
    291291    psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row
    292     if (!effGenerate(&xFake, &yFake, readout, psf, magOffsets,
    293                      numSources, magLim, radius, minFlux)) {
     292    pmReadout *fakeRO = effGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius, minFlux);
     293    if (!fakeRO) {
    294294        psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources");
    295295        psFree(xFake);
    … …  
    416416    psFree(significance);
    417417
     418    // psphotFitSourcesLinearReadout subtracts the model fits
    418419    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
    419420        psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
    … …  
    427428    psf->ApTrend = NULL;
    428429
     430    // measure the magnitudes and fluxes for the sources
    429431    if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf)) {
    430432        psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources.");
    … …  
    433435        psf->ApTrend = apTrend; // Casting away const!
    434436        return false;
     437    }
     438
     439    // replace the subtracted model fits
     440    for (int i = 0; i < fakeSourcesAll->n; i++) {
     441        pmSource *source = fakeSourcesAll->data[i];
     442
     443        // replace other sources?
     444        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
     445       
     446        pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    435447    }
    436448    psFree(fakeSourcesAll);
    … …  
    515527    psFree(fakeSources);
    516528
     529    // subtract the faked sources from the original image
     530    psBinaryOp(readout->image, readout->image, "-", fakeRO->image);
     531    psFree(fakeRO);
     532
    517533    pmDetEff *de = pmDetEffAlloc(magLim, numSources, numBins); // Detection efficiency
    518534    de->magOffsets = psVectorCopy(NULL, magOffsets, PS_TYPE_F32);

  • branches/eam_branches/ipp-20110710/psphot/src/psphotEllipticalContour.c

    r29004 r32023  
    127127psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord) {
    128128
    129     static int pass = 0;
    130 
    131129    psF32 *par = params->data.F32;
    132130
    … …  
    145143
    146144    // value is X
    147     // if (coord->data.F32[1] == 0) {
    148     if (pass == 0) {
    149         pass = 1;
    150 
     145    if (coord->data.F32[1] < 0.5) {
    151146        float value = par[PAR_RMIN]*cs_alpha*r;
    152147
    … …  
    161156
    162157    // value is Y
    163     // if (coord->data.F32[1] == 1) {
    164     if (pass == 1) {
    165         pass = 0;
    166 
     158    if (coord->data.F32[1] > 0.5) {
    167159        float value = par[PAR_RMIN]*sn_alpha*r;
    168160

  • branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceAnalysis.c

    r31154 r32023  
    22
    33// measure the elliptical radial profile and use this to measure the petrosian parameters for the sources
    4 // XXX this function needs to be threaded
    54
    65// for now, let's store the detections on the readout->analysis for each readout
    … …  
    5958    }
    6059
    61     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    62     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    63     assert (maskVal);
     60    // determine the number of allowed threads
     61    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     62    if (!status) {
     63        nThreads = 0;
     64    }
    6465
    6566    // get the sky noise from the background analysis; if this is missing, get the user-supplied value
    … …  
    7071    }
    7172
     73    // source analysis is done in S/N order (brightest first)
     74    sources = psArraySort (sources, pmSourceSortByFlux);
     75
     76    // option to limit analysis to a specific region
     77    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
     78    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
     79    *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
     80    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
     81
     82    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
     83    int Cx = 1, Cy = 1;
     84    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
     85
     86    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
     87
     88    for (int i = 0; i < cellGroups->n; i++) {
     89
     90        psArray *cells = cellGroups->data[i];
     91
     92        for (int j = 0; j < cells->n; j++) {
     93
     94            // allocate a job -- if threads are not defined, this just runs the job
     95            psThreadJob *job = psThreadJobAlloc ("PSPHOT_EXTENDED_ANALYSIS");
     96
     97            psArrayAdd(job->args, 1, readout);
     98            psArrayAdd(job->args, 1, cells->data[j]); // sources
     99            psArrayAdd(job->args, 1, AnalysisRegion);
     100            psArrayAdd(job->args, 1, recipe);
     101
     102            PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32);
     103
     104            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
     105            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
     106            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
     107
     108// set this to 0 to run without threading
     109# if (1)           
     110            if (!psThreadJobAddPending(job)) {
     111                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     112                psFree(AnalysisRegion);
     113                return false;
     114            }
     115# else
     116            if (!psphotExtendedSourceAnalysis_Threaded(job)) {
     117                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     118                psFree(AnalysisRegion);
     119                return false;
     120            }
     121            psScalar *scalar = NULL;
     122            scalar = job->args->data[5];
     123            Next += scalar->data.S32;
     124            scalar = job->args->data[6];
     125            Npetro += scalar->data.S32;
     126            scalar = job->args->data[7];
     127            Nannuli += scalar->data.S32;
     128            psFree(job);
     129# endif
     130        }
     131
     132        // wait for the threads to finish and manage results
     133        if (!psThreadPoolWait (false)) {
     134            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     135            psFree(AnalysisRegion);
     136            return false;
     137        }
     138
     139        // we have only supplied one type of job, so we can assume the types here
     140        psThreadJob *job = NULL;
     141        while ((job = psThreadJobGetDone()) != NULL) {
     142            if (job->args->n < 1) {
     143                fprintf (stderr, "error with job\n");
     144            } else {
     145                psScalar *scalar = NULL;
     146                scalar = job->args->data[5];
     147                Next += scalar->data.S32;
     148                scalar = job->args->data[6];
     149                Npetro += scalar->data.S32;
     150                scalar = job->args->data[7];
     151                Nannuli += scalar->data.S32;
     152            }
     153            psFree(job);
     154        }
     155    }
     156    psFree (cellGroups);
     157    psFree(AnalysisRegion);
     158
     159    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     160    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
     161    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
     162
     163    psphotVisualShowResidualImage (readout, false);
     164
     165    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     166    if (doPetrosian) {
     167        psphotVisualShowPetrosians (sources);
     168    }
     169
     170    return true;
     171}
     172
     173bool psphotExtendedSourceAnalysis_Threaded (psThreadJob *job) {
     174
     175    bool status;
     176
     177    int Next = 0;
     178    int Npetro = 0;
     179    int Nannuli = 0;
     180
     181    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
     182    pmReadout *readout      = job->args->data[0];
     183    psArray *sources        = job->args->data[1];
     184    psRegion *region        = job->args->data[2];
     185    psMetadata *recipe      = job->args->data[3];
     186
     187    float skynoise          = PS_SCALAR_VALUE(job->args->data[4],F32);
     188
    72189    // S/N limit to perform full non-linear fits
    73190    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
    … …  
    78195    bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
    79196
    80     // source analysis is done in S/N order (brightest first)
    81     sources = psArraySort (sources, pmSourceSortByFlux);
    82 
    83     // option to limit analysis to a specific region
    84     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    85     psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    86     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
     197    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     198    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     199    assert (maskVal);
    87200
    88201    // choose the sources of interest
    … …  
    108221
    109222        // limit selection by analysis region
    110         if (source->peak->x < AnalysisRegion.x0) continue;
    111         if (source->peak->y < AnalysisRegion.y0) continue;
    112         if (source->peak->x > AnalysisRegion.x1) continue;
    113         if (source->peak->y > AnalysisRegion.y1) continue;
     223        if (source->peak->x < region->x0) continue;
     224        if (source->peak->y < region->y0) continue;
     225        if (source->peak->x > region->x1) continue;
     226        if (source->peak->y > region->y1) continue;
    114227
    115228        // replace object in image
    … …  
    159272    }
    160273
    161     psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
    162     psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
    163     psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
    164 
    165     psphotVisualShowResidualImage (readout, false);
    166 
    167     if (doPetrosian) {
    168         psphotVisualShowPetrosians (sources);
    169     }
    170 
    171     // fprintf (stderr, "xsrc : tried %ld objects\n", sources->n);
     274    psScalar *scalar = NULL;
     275
     276    // change the value of a scalar on the array (wrap this and put it in psArray.h)
     277    scalar = job->args->data[5];
     278    scalar->data.S32 = Next;
     279
     280    scalar = job->args->data[6];
     281    scalar->data.S32 = Npetro;
     282
     283    scalar = job->args->data[7];
     284    scalar->data.S32 = Nannuli;
    172285
    173286    return true;

  • branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c

    r31990 r32023  
    278278    psMetadata *models      = job->args->data[2];
    279279    psRegion *region        = job->args->data[3];
    280     int psfSize             = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     280    int psfSize             = PS_SCALAR_VALUE(job->args->data[4],S32);
    281281    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
    282282    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);

  • branches/eam_branches/ipp-20110710/psphot/src/psphotRadialApertures.c

    r31452 r32023  
    6363    }
    6464
     65    // determine the number of allowed threads
     66    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     67    if (!status) {
     68        nThreads = 0;
     69    }
     70
     71    // source analysis is done in S/N order (brightest first)
     72    // XXX are we getting the objects out of order? does it matter?
     73    sources = psArraySort (sources, pmSourceSortByFlux);
     74
     75    // option to limit analysis to a specific region
     76    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
     77    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
     78    *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
     79    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
     80
     81    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
     82    int Cx = 1, Cy = 1;
     83    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
     84
     85    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
     86
     87    for (int i = 0; i < cellGroups->n; i++) {
     88
     89        psArray *cells = cellGroups->data[i];
     90
     91        for (int j = 0; j < cells->n; j++) {
     92
     93            // allocate a job -- if threads are not defined, this just runs the job
     94            psThreadJob *job = psThreadJobAlloc ("PSPHOT_RADIAL_APERTURES");
     95
     96            psArrayAdd(job->args, 1, readout);
     97            psArrayAdd(job->args, 1, cells->data[j]); // sources
     98            psArrayAdd(job->args, 1, AnalysisRegion);
     99            psArrayAdd(job->args, 1, recipe);
     100
     101            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
     102
     103// set this to 0 to run without threading
     104# if (1)           
     105            if (!psThreadJobAddPending(job)) {
     106                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     107                psFree(AnalysisRegion);
     108                return false;
     109            }
     110# else
     111            if (!psphotRadialApertures_Threaded(job)) {
     112                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     113                psFree(AnalysisRegion);
     114                return false;
     115            }
     116            psScalar *scalar = NULL;
     117            scalar = job->args->data[4];
     118            Nradial += scalar->data.S32;
     119            psFree(job);
     120# endif
     121        }
     122
     123        // wait for the threads to finish and manage results
     124        if (!psThreadPoolWait (false)) {
     125            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     126            psFree(AnalysisRegion);
     127            return false;
     128        }
     129
     130        // we have only supplied one type of job, so we can assume the types here
     131        psThreadJob *job = NULL;
     132        while ((job = psThreadJobGetDone()) != NULL) {
     133            if (job->args->n < 1) {
     134                fprintf (stderr, "error with job\n");
     135            } else {
     136                psScalar *scalar = NULL;
     137                scalar = job->args->data[4];
     138                Nradial += scalar->data.S32;
     139            }
     140            psFree(job);
     141        }
     142    }
     143    psFree (cellGroups);
     144    psFree(AnalysisRegion);
     145
     146    psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
     147    return true;
     148}
     149 
     150bool psphotRadialApertures_Threaded (psThreadJob *job) {
     151
     152    bool status;
     153    int Nradial = 0;
     154
     155    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
     156    pmReadout *readout      = job->args->data[0];
     157    psArray *sources        = job->args->data[1];
     158    psRegion *region        = job->args->data[2];
     159    psMetadata *recipe      = job->args->data[3];
     160
    65161    // radMax stores the upper bounds of the annuli
    66162    // XXX keep the same name here as for the petrosian / elliptical apertures?
    … …  
    68164    psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
    69165    psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
    70     float outerRadius = radMax->data.F32[radMax->n - 1];
    71166
    72167    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    77172    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
    78173
    79     // source analysis is done in S/N order (brightest first)
    80     // XXX are we getting the objects out of order? does it matter?
    81     sources = psArraySort (sources, pmSourceSortByFlux);
    82 
    83     // option to limit analysis to a specific region
    84     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    85     psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    86     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
     174    float outerRadius = radMax->data.F32[radMax->n - 1];
    87175
    88176    // choose the sources of interest
    … …  
    104192
    105193        // limit selection by analysis region
    106         if (source->peak->x < AnalysisRegion.x0) continue;
    107         if (source->peak->y < AnalysisRegion.y0) continue;
    108         if (source->peak->x > AnalysisRegion.x1) continue;
    109         if (source->peak->y > AnalysisRegion.y1) continue;
     194        if (source->peak->x < region->x0) continue;
     195        if (source->peak->y < region->y0) continue;
     196        if (source->peak->x > region->x1) continue;
     197        if (source->peak->y > region->y1) continue;
    110198
    111199        // allocate pmSourceExtendedParameters, if not already defined
    … …  
    145233        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    146234    }
    147 
    148     psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
     235    psScalar *scalar = job->args->data[4];
     236    scalar->data.S32 = Nradial;
     237
    149238    return true;
    150239}

  • branches/eam_branches/ipp-20110710/psphot/src/psphotSetThreads.c

    r31452 r32023  
    4040    psFree(task);
    4141
     42    task = psThreadTaskAlloc("PSPHOT_EXTENDED_ANALYSIS", 8);
     43    task->function = &psphotExtendedSourceAnalysis_Threaded;
     44    psThreadTaskAdd(task);
     45    psFree(task);
     46
     47    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 5);
     48    task->function = &psphotRadialApertures_Threaded;
     49    psThreadTaskAdd(task);
     50    psFree(task);
     51
    4252    return true;
    4353}

  • branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c

    r31990 r32023  
    637637
    638638        pmSourceFitModel (source, model, &options, maskVal);
    639         fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
     639        // fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    640640
    641641        chiSquare[i] = model->chisqNorm;
Note: See TracChangeset for help on using the changeset viewer.