Index: /trunk/ippconfig/recipes/psphot.config =================================================================== --- /trunk/ippconfig/recipes/psphot.config (revision 37465) +++ /trunk/ippconfig/recipes/psphot.config (revision 37466) @@ -233,7 +233,8 @@ # Extended source fit parameters # these models are used for high-quality shape analysis after all sources have been removed -EXTENDED_SOURCE_FITS BOOL FALSE # perform the detailed extended sourc fits? -EXTENDED_SOURCE_FITS_INTERACTIVE BOOL FALSE # perform the detailed extended sourc fits? +EXTENDED_SOURCE_FITS BOOL FALSE # perform the detailed extended source fits? +EXTENDED_SOURCE_FITS_INTERACTIVE BOOL FALSE # perform the detailed extended source fits? EXTENDED_SOURCE_FITS_POISSON BOOL FALSE # use poisson errors for extended source models +EXTENDED_SOURCE_GALAXY_PARAMS BOOL FALSE # measure galaxy parameters EXTENDED_SOURCE_MODELS_SELECTION STR ALL @@ -479,4 +480,5 @@ STACKPHOT METADATA EXTENDED_SOURCE_FITS BOOL TRUE # perform any of the aperture-like measurements? + EXTENDED_SOURCE_GALAXY_PARAMS BOOL TRUE # perform any of the aperture-like measurements? EXTENDED_SOURCE_SN_LIM F32 10.0 EXTENDED_SOURCE_PETROSIAN BOOL TRUE @@ -494,5 +496,5 @@ RADIAL_APERTURES BOOL T # calculate flux in circular radial apertures? RADIAL_APERTURES_SN_LIM F32 0.0 # S/N limit for radial aperture calculation - OUTPUT.FORMAT STR PS1_SV3 + OUTPUT.FORMAT STR PS1_SV4 # subtract radial profiles for extremely bright stars? Index: /trunk/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- /trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 37465) +++ /trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 37466) @@ -226,4 +226,10 @@ pars->petrosianR50 = NAN; pars->petrosianR50Err = NAN; + pars->ghalfLightRadius = NAN; + pars->gRT = NAN; + pars->gRA = NAN; + pars->gS2 = NAN; + pars->gA = NAN; + pars->gbumpy = NAN; return pars; } Index: /trunk/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- /trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 37465) +++ /trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 37466) @@ -65,4 +65,10 @@ float petrosianR50Err; float petrosianFill; + float ghalfLightRadius; + float gRT; // total residual (from eliptically symmetric model) Simard 2002 + float gRA; // assymetric residual (from eliptically symmetric model) Simard 2002 + float gS2; // (un) smoothness (Simard, 2011 Cheng 2011) + float gA; // assymetry index (Gyory & Bell 2010) + float gbumpy; // Blakeslee bumpiness } pmSourceExtendedPars; Index: /trunk/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 37465) +++ /trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 37466) @@ -639,4 +639,10 @@ psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "HALF_LIGHT_RADIUS", 0, "Half Light Radius", extpars->ghalfLightRadius); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "G_RT", 0, "total residual Sersic Model", extpars->gRT); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "G_RA", 0, "assymetric residual Sersic Model", extpars->gRA); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "G_S2", 0, "(un)smoothness paramter residual Sersic Model", extpars->gS2); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "G_A", 0, "assymetry index", extpars->gA); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "G_BUMPY", 0, "bumpiness", extpars->gbumpy); } else { psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); @@ -649,4 +655,10 @@ psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", NAN); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "HALF_LIGHT_RADIUS", 0, "Half Light Radius", NAN); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "G_RT", 0, "total residual Sersic Model", NAN); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "G_RA", 0, "assymetric residual Sersic Model", NAN); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "G_S2", 0, "(un)smoothness paramter residual Sersic Model", NAN); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "G_A", 0, "assymetry index", NAN); + @>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "G_BUMPY", 0, "bumpiness", NAN); } } Index: /trunk/psphot/src/Makefile.am =================================================================== --- /trunk/psphot/src/Makefile.am (revision 37465) +++ /trunk/psphot/src/Makefile.am (revision 37466) @@ -260,5 +260,6 @@ psphotSetNFrames.c \ psphotSourceMemory.c \ - psphotChipParams.c + psphotChipParams.c \ + psphotGalaxyParams.c # not currently used Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 37465) +++ /trunk/psphot/src/psphot.h (revision 37466) @@ -137,4 +137,8 @@ bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); bool psphotExtendedSourceFits_Threaded (psThreadJob *job); + +bool psphotGalaxyParams (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotGalaxyParamsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotGalaxyParams_Threaded (psThreadJob *job); bool psphotApResid (pmConfig *config, const pmFPAview *view, const char *filerule); Index: /trunk/psphot/src/psphotGalaxyParams.c =================================================================== --- /trunk/psphot/src/psphotGalaxyParams.c (revision 37466) +++ /trunk/psphot/src/psphotGalaxyParams.c (revision 37466) @@ -0,0 +1,494 @@ +# include "psphotInternal.h" + +// measure some properties on the exended objects + +bool psphotGalaxyParams (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Galaxy Paramters ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // add this + bool doGalaxyParams = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_GALAXY_PARAMS"); + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doExtModels = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS"); + + // All of these must be true to measure galaxy parameters + if (!doGalaxyParams || !doPetrosian || !doAnnuli || !doExtModels) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping Galaxy parameter source measurements\n"); + return true; + } + + int num = psphotFileruleCount(config, filerule); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotGalaxyParamsReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on galaxy paramter measurements for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +#ifdef notdef +/*** for the moment, this test code : it is not thread safe ***/ +int Nall = 0; +int Nskip1 = 0; +int Nskip2 = 0; +int Nskip3 = 0; +int Nskip4 = 0; +int Nskip5 = 0; +int Nskip6 = 0; + +# define SKIP(VALUE) { VALUE++; continue; } + +#endif //notdef + +#define SKIP(VALUE) {continue;} + +// aperture-like measurements for extended sources +bool psphotGalaxyParamsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { + + bool status; + int Next = 0; + int Npetro = 0; + int Nannuli = 0; + + psTimerStart ("psphot.galaxy"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + return true; + } + + // 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 + float skynoise = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV"); + if (!status) { + skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + psWarning ("failed to get sky noise level from background analysis; defaulting to user supplied value of %f\n", skynoise); + } + + // 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_GALAXY_PARAMS"); + + 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 (0) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } +# else + if (!psphotGalaxyParams_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, true)) { + 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_WARN, "galaxy parameter analysis: %f sec for %d objects\n", psTimerMark ("psphot.galaxy"), Next); + psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); + psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); + +// # if (PS_TRACE_ON) +# if (0) + fprintf (stderr, "ext analysis skipped @ 1 : %d\n", Nskip1); + fprintf (stderr, "ext analysis skipped @ 2 : %d\n", Nskip2); + fprintf (stderr, "ext analysis skipped @ 3 : %d\n", Nskip3); + fprintf (stderr, "ext analysis skipped @ 4 : %d\n", Nskip4); + fprintf (stderr, "ext analysis skipped @ 5 : %d\n", Nskip5); + fprintf (stderr, "ext analysis skipped @ 6 : %d\n", Nskip6); +#endif + + // psphotRadialProfileShowSkips (); + + // psphotVisualShowResidualImage (readout, false); + + // psphotVisualShowPetrosians (sources); + + return true; +} + +bool psphotGalaxyParams_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); + + bool status; +// bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + + // 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); + + pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + + psVector *bgPixels = NULL; + psVector *bgPixelsFlip = NULL; +#ifdef MAKE_MARK_IMAGE + psImage *markImage = NULL; +#endif + + #define MARK_SOURCE 0 + #define MARK_MASK 1 + #define MARK_BG 2 + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + + // choose the sources of interest + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // Nall ++; + + // rules for measuring petrosian parameters for specific objects are set in + // psphotChooseAnalysisOptions.c + if (!(source->tmpFlags & PM_SOURCE_TMPF_PETRO)) SKIP (Nskip1); + + // limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?) + if (source->peak->x < region->x0) SKIP (Nskip2); + if (source->peak->y < region->y0) SKIP (Nskip3); + if (source->peak->x > region->x1) SKIP (Nskip4); + if (source->peak->y > region->y1) SKIP (Nskip5); + + // Should we skip objects that do not have an extended model? + if (!source->extpars) SKIP (Nskip6); + if (!source->modelFits) SKIP (Nskip7); + + pmModel *sersicModel = NULL; + for (int j = 0; j < source->modelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + if (model->type == sersicModelType) { + sersicModel = model; + break; + } + } + if (!sersicModel) SKIP (Nskip8); + Next ++; + + // if sersicModel is not modelEXT make it so + pmSourceMode saveMode = source->mode; + pmSourceType saveType = source->type; + bool isExtended = saveType == PM_SOURCE_TYPE_EXTENDED; + pmModel *extModel = source->modelEXT; + + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + if (!isExtended || extModel != sersicModel) { + source->modelEXT = psMemIncrRefCounter (sersicModel); + source->type = PM_SOURCE_TYPE_EXTENDED; + source->mode |= PM_SOURCE_MODE_EXTMODEL; + // do we need this? + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; + } + + // force source image to be a bit larger... + float radius = source->peak->xf - source->pixels->col0; + radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius); + + // Do the work here + + psF32 *PAR = sersicModel->params->data.F32; + psF32 Xo = PAR[PM_PAR_XPOS]; + psF32 Yo = PAR[PM_PAR_YPOS]; + + // For now: set the half light radius to the major axis of the model + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, sersicModel->class->useReff); + source->extpars->ghalfLightRadius = axes.major; + psF32 R = 2 * axes.major; + psF32 R2 = PS_SQR(R); + + // Create the markImage and accumlate the Flux + + // center of mass in subimage. Note: the calculation below uses pixel index, so we correct + // xCM, yCM from pixel coords to pixel index here. + psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage + psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage + + psF32 **vPix = source->pixels->data.F32; + psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + + bgPixels = psVectorRecycle(bgPixels, source->pixels->numCols*source->pixels->numRows, PS_TYPE_F32); + bgPixels->n = 0; + bgPixelsFlip = psVectorRecycle(bgPixelsFlip, source->pixels->numCols*source->pixels->numRows, PS_TYPE_F32); + bgPixelsFlip->n = 0; + +#ifdef MAKE_MARK_IMAGE + // The markImage may no longer be needed + markImage = psImageRecycle(markImage, source->pixels->numCols, source->pixels->numRows, PS_TYPE_U8); + + // assume all pixels are masked to start + psImageInit(markImage, MARK_MASK); + psU8 **vMark = markImage->data.U8; + #define SET_MARK(_val) vMark[row][col] = _val +#else + #define SET_MARK(_val) +#endif + + // XXX use recipe value instead of hard coded 1.5 + psF32 threshold = skynoise * 1.5; + + psF32 sumI = 0; + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + psF32 yDiff = row - yCM; + // y coordinate of mirror pixel + int yFlip = yCM - yDiff; + if (yFlip < 0) continue; + if (yFlip >= source->pixels->numRows) continue; + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + // check mask and value for this pixel + if (vMsk && (vMsk[row][col] & maskVal)) continue; + psF32 pixel = vPix[row][col]; + if (isnan(pixel)) continue; + // valid pixel + + psF32 xDiff = col - xCM; + + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 <= R2) { + sumI += pixel; + } + + + if (fabs(pixel) >= threshold) { + SET_MARK(MARK_SOURCE); + continue; + } + // we have a background pixel. + // Check whether it's mirror is a background pixel as well. + + // x coordinate of mirror pixel + int xFlip = xCM - xDiff; + if (xFlip < 0) continue; + if (xFlip >= source->pixels->numCols) continue; + // check mask and value for mirror pixel + if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue; + psF32 mirrorPixel = vPix[yFlip][xFlip]; + if (isnan(mirrorPixel)) continue; + if (fabs(mirrorPixel) < threshold) { + // This is a good background pixel with a good mirror Save the values + SET_MARK(MARK_BG); + psVectorAppend(bgPixels, pixel); + psVectorAppend(bgPixelsFlip, mirrorPixel); + } + } + } + + // subtract model to create the residual image + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { + // HMM something has gone wrong + goto restoreModel; + } + + // don't know if these could change in the background subtraction + vPix = source->pixels->data.F32; + vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + // Accumulate Rt and Ra + psF32 sumRt = 0; + psF32 sumRa = 0; + int numPixels = 0; + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + psF32 yDiff = row - yCM; + // y coordinate of mirror pixel + int yFlip = yCM - yDiff; + if (yFlip < 0) continue; + if (yFlip >= source->pixels->numRows) continue; + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + // check mask and value for this pixel + if (vMsk && (vMsk[row][col] & maskVal)) continue; + psF32 pixel = vPix[row][col]; + if (isnan(pixel)) continue; + + psF32 xDiff = col - xCM; + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > R2) continue; + + // x coordinate of mirror pixel + int xFlip = xCM - xDiff; + if (xFlip < 0) continue; + if (xFlip >= source->pixels->numCols) continue; + // check mask and value for mirror pixel + if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue; + psF32 mirrorPixel = vPix[yFlip][xFlip]; + if (isnan(mirrorPixel)) continue; + + sumRt += fabs(pixel + mirrorPixel); + sumRa += fabs(pixel - mirrorPixel); + ++numPixels; + } + } + // and calculate Bt and Ba + psF32 sumBt = 0; + psF32 sumBa = 0; + psF32 *vB = bgPixels->data.F32; + psF32 *vBFlip = bgPixelsFlip->data.F32; + long length = bgPixels->n; + // Randomly select numPixels background pixels + for (int n = 0; n < numPixels; n++) { + double rand = psRandomUniform(rng); + long index = floor(rand * length); + if (index >= length) continue; + sumBt += fabs(vB[index] + vBFlip[index]); + sumBa += fabs(vB[index] - vBFlip[index]); + } + + // finally save the results + source->extpars->gRT = 0.5 * (sumRt - sumBt)/sumI; + source->extpars->gRA = 0.5 * (sumRa - sumBa)/sumI; + source->extpars->gS2 = source->extpars->gRT + source->extpars->gRA; + + + // Now put things back the way that they were. +restoreModel: + + if ((extModel != sersicModel) || !isExtended) { + // add the sersic model back in + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + // replace the extended model + if (extModel != sersicModel) { + psMemDecrRefCounter (sersicModel); + source->modelEXT = extModel; + } + source->mode = saveMode; + source->type = saveType; + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + + } + +#ifdef MAKE_MARK_IMAGE + psFree(markImage); +#endif + psFree(bgPixels) + psFree(bgPixelsFlip) + psFree(rng); + + 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: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 37465) +++ /trunk/psphot/src/psphotReadout.c (revision 37466) @@ -325,4 +325,6 @@ psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources) psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources) + // measure some parameters for galaxy science + psphotGalaxyParams (config, view, filerule); psphotRadialApertures(config, view, filerule, 0); Index: /trunk/psphot/src/psphotSetThreads.c =================================================================== --- /trunk/psphot/src/psphotSetThreads.c (revision 37465) +++ /trunk/psphot/src/psphotSetThreads.c (revision 37466) @@ -92,4 +92,9 @@ psFree(task); + task = psThreadTaskAlloc("PSPHOT_GALAXY_PARAMS", 8); + task->function = &psphotGalaxyParams_Threaded; + psThreadTaskAdd(task); + psFree(task); + return true; } Index: /trunk/psphot/src/psphotStackReadout.c =================================================================== --- /trunk/psphot/src/psphotStackReadout.c (revision 37465) +++ /trunk/psphot/src/psphotStackReadout.c (revision 37466) @@ -342,4 +342,6 @@ psphotExtendedSourceFits (config, view, STACK_RAW); // pass 1 (detections->allSources) + // measure some parameters for galxy science + psphotGalaxyParams (config, view, STACK_RAW); // pass 1 (detections->allSources) // create source children for the OUT filerule (for radial aperture photometry and output) // NOTE: The new source children have image arrays pointing to the readout associated with