Index: /branches/eam_branches/ipp-20110710/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphot.h (revision 32022) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphot.h (revision 32023) @@ -117,4 +117,5 @@ bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotExtendedSourceAnalysis_Threaded (psThreadJob *job); bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule); @@ -445,4 +446,5 @@ bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotRadialApertures_Threaded (psThreadJob *job); bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry); Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotEfficiency.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotEfficiency.c (revision 32022) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotEfficiency.c (revision 32023) @@ -99,5 +99,5 @@ /// Generate a fake image and add it in to the existing readout -static bool effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources +static pmReadout *effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources const pmReadout *ro, // Readout of interest const pmPSF *psf, // Point-spread function @@ -152,5 +152,5 @@ psFree(xAll); psFree(yAll); - return false; + return NULL; } psFree(magAll); @@ -161,7 +161,7 @@ psBinaryOp(ro->image, ro->image, "+", fakeRO->image); - psFree(fakeRO); - - return true; + + // return the readout so we can subtract it later + return fakeRO; } @@ -290,6 +290,6 @@ psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row - if (!effGenerate(&xFake, &yFake, readout, psf, magOffsets, - numSources, magLim, radius, minFlux)) { + pmReadout *fakeRO = effGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius, minFlux); + if (!fakeRO) { psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources"); psFree(xFake); @@ -416,4 +416,5 @@ psFree(significance); + // psphotFitSourcesLinearReadout subtracts the model fits if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); @@ -427,4 +428,5 @@ psf->ApTrend = NULL; + // measure the magnitudes and fluxes for the sources if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources."); @@ -433,4 +435,14 @@ psf->ApTrend = apTrend; // Casting away const! return false; + } + + // replace the subtracted model fits + for (int i = 0; i < fakeSourcesAll->n; i++) { + pmSource *source = fakeSourcesAll->data[i]; + + // replace other sources? + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; + + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } psFree(fakeSourcesAll); @@ -515,4 +527,8 @@ psFree(fakeSources); + // subtract the faked sources from the original image + psBinaryOp(readout->image, readout->image, "-", fakeRO->image); + psFree(fakeRO); + pmDetEff *de = pmDetEffAlloc(magLim, numSources, numBins); // Detection efficiency de->magOffsets = psVectorCopy(NULL, magOffsets, PS_TYPE_F32); Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotEllipticalContour.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotEllipticalContour.c (revision 32022) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotEllipticalContour.c (revision 32023) @@ -127,6 +127,4 @@ psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord) { - static int pass = 0; - psF32 *par = params->data.F32; @@ -145,8 +143,5 @@ // value is X - // if (coord->data.F32[1] == 0) { - if (pass == 0) { - pass = 1; - + if (coord->data.F32[1] < 0.5) { float value = par[PAR_RMIN]*cs_alpha*r; @@ -161,8 +156,5 @@ // value is Y - // if (coord->data.F32[1] == 1) { - if (pass == 1) { - pass = 0; - + if (coord->data.F32[1] > 0.5) { float value = par[PAR_RMIN]*sn_alpha*r; Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceAnalysis.c (revision 32022) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceAnalysis.c (revision 32023) @@ -2,5 +2,4 @@ // measure the elliptical radial profile and use this to measure the petrosian parameters for the sources -// XXX this function needs to be threaded // for now, let's store the detections on the readout->analysis for each readout @@ -59,7 +58,9 @@ } - // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) - psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels - assert (maskVal); + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } // get the sky noise from the background analysis; if this is missing, get the user-supplied value @@ -70,4 +71,120 @@ } + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortByFlux); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0); + *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined"); + + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_EXTENDED_ANALYSIS"); + + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, AnalysisRegion); + psArrayAdd(job->args, 1, recipe); + + PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32); + + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli + +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } +# else + if (!psphotExtendedSourceAnalysis_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } + psScalar *scalar = NULL; + scalar = job->args->data[5]; + Next += scalar->data.S32; + scalar = job->args->data[6]; + Npetro += scalar->data.S32; + scalar = job->args->data[7]; + Nannuli += scalar->data.S32; + psFree(job); +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[5]; + Next += scalar->data.S32; + scalar = job->args->data[6]; + Npetro += scalar->data.S32; + scalar = job->args->data[7]; + Nannuli += scalar->data.S32; + } + psFree(job); + } + } + psFree (cellGroups); + psFree(AnalysisRegion); + + psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); + psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); + psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); + + psphotVisualShowResidualImage (readout, false); + + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + if (doPetrosian) { + psphotVisualShowPetrosians (sources); + } + + return true; +} + +bool psphotExtendedSourceAnalysis_Threaded (psThreadJob *job) { + + bool status; + + int Next = 0; + int Npetro = 0; + int Nannuli = 0; + + // arguments: readout, sources, models, region, psfSize, maskVal, markVal + pmReadout *readout = job->args->data[0]; + psArray *sources = job->args->data[1]; + psRegion *region = job->args->data[2]; + psMetadata *recipe = job->args->data[3]; + + float skynoise = PS_SCALAR_VALUE(job->args->data[4],F32); + // S/N limit to perform full non-linear fits float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); @@ -78,11 +195,7 @@ bool doPetroStars = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS"); - // source analysis is done in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortByFlux); - - // option to limit analysis to a specific region - char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); - psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); - if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); // choose the sources of interest @@ -108,8 +221,8 @@ // limit selection by analysis region - if (source->peak->x < AnalysisRegion.x0) continue; - if (source->peak->y < AnalysisRegion.y0) continue; - if (source->peak->x > AnalysisRegion.x1) continue; - if (source->peak->y > AnalysisRegion.y1) continue; + if (source->peak->x < region->x0) continue; + if (source->peak->y < region->y0) continue; + if (source->peak->x > region->x1) continue; + if (source->peak->y > region->y1) continue; // replace object in image @@ -159,15 +272,15 @@ } - psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); - psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); - psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); - - psphotVisualShowResidualImage (readout, false); - - if (doPetrosian) { - psphotVisualShowPetrosians (sources); - } - - // fprintf (stderr, "xsrc : tried %ld objects\n", sources->n); + psScalar *scalar = NULL; + + // change the value of a scalar on the array (wrap this and put it in psArray.h) + scalar = job->args->data[5]; + scalar->data.S32 = Next; + + scalar = job->args->data[6]; + scalar->data.S32 = Npetro; + + scalar = job->args->data[7]; + scalar->data.S32 = Nannuli; return true; Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c (revision 32022) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c (revision 32023) @@ -278,5 +278,5 @@ psMetadata *models = job->args->data[2]; psRegion *region = job->args->data[3]; - int psfSize = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); + int psfSize = PS_SCALAR_VALUE(job->args->data[4],S32); psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotRadialApertures.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotRadialApertures.c (revision 32022) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotRadialApertures.c (revision 32023) @@ -63,4 +63,100 @@ } + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } + + // source analysis is done in S/N order (brightest first) + // XXX are we getting the objects out of order? does it matter? + sources = psArraySort (sources, pmSourceSortByFlux); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0); + *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined"); + + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_RADIAL_APERTURES"); + + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, AnalysisRegion); + psArrayAdd(job->args, 1, recipe); + + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial + +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } +# else + if (!psphotRadialApertures_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } + psScalar *scalar = NULL; + scalar = job->args->data[4]; + Nradial += scalar->data.S32; + psFree(job); +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[4]; + Nradial += scalar->data.S32; + } + psFree(job); + } + } + psFree (cellGroups); + psFree(AnalysisRegion); + + psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial); + return true; +} + +bool psphotRadialApertures_Threaded (psThreadJob *job) { + + bool status; + int Nradial = 0; + + // arguments: readout, sources, models, region, psfSize, maskVal, markVal + pmReadout *readout = job->args->data[0]; + psArray *sources = job->args->data[1]; + psRegion *region = job->args->data[2]; + psMetadata *recipe = job->args->data[3]; + // radMax stores the upper bounds of the annuli // XXX keep the same name here as for the petrosian / elliptical apertures? @@ -68,5 +164,4 @@ psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe"); psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define"); - float outerRadius = radMax->data.F32[radMax->n - 1]; // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -77,12 +172,5 @@ float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM"); - // source analysis is done in S/N order (brightest first) - // XXX are we getting the objects out of order? does it matter? - sources = psArraySort (sources, pmSourceSortByFlux); - - // option to limit analysis to a specific region - char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); - psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); - if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + float outerRadius = radMax->data.F32[radMax->n - 1]; // choose the sources of interest @@ -104,8 +192,8 @@ // limit selection by analysis region - if (source->peak->x < AnalysisRegion.x0) continue; - if (source->peak->y < AnalysisRegion.y0) continue; - if (source->peak->x > AnalysisRegion.x1) continue; - if (source->peak->y > AnalysisRegion.y1) continue; + if (source->peak->x < region->x0) continue; + if (source->peak->y < region->y0) continue; + if (source->peak->x > region->x1) continue; + if (source->peak->y > region->y1) continue; // allocate pmSourceExtendedParameters, if not already defined @@ -145,6 +233,7 @@ pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); } - - psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial); + psScalar *scalar = job->args->data[4]; + scalar->data.S32 = Nradial; + return true; } Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotSetThreads.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotSetThreads.c (revision 32022) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotSetThreads.c (revision 32023) @@ -40,4 +40,14 @@ psFree(task); + task = psThreadTaskAlloc("PSPHOT_EXTENDED_ANALYSIS", 8); + task->function = &psphotExtendedSourceAnalysis_Threaded; + psThreadTaskAdd(task); + psFree(task); + + task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 5); + task->function = &psphotRadialApertures_Threaded; + psThreadTaskAdd(task); + psFree(task); + return true; } Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c (revision 32022) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c (revision 32023) @@ -637,5 +637,5 @@ pmSourceFitModel (source, model, &options, maskVal); - fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix); + // fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix); chiSquare[i] = model->chisqNorm;