Index: trunk/psphot/src/psphotGalaxyParams.c =================================================================== --- trunk/psphot/src/psphotGalaxyParams.c (revision 37625) +++ trunk/psphot/src/psphotGalaxyParams.c (revision 37626) @@ -1,3 +1,5 @@ # include "psphotInternal.h" + +#define USE_FOOTPRINTS 1 // measure some properties on the exended objects @@ -95,4 +97,18 @@ } + float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolutio + if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) { + fitNsigmaConv = 5.0; + } + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + + psImage *footprintImage = NULL; +#ifdef USE_FOOTPRINTS + footprintImage = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_S32); + psImageInit(footprintImage, 0); + pmSetFootprintArrayIDsForImage(footprintImage, detections->footprints, true); +#endif + + // source analysis is done in S/N order (brightest first) sources = psArraySort (sources, pmSourceSortByFlux); @@ -125,8 +141,12 @@ PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, fitNsigmaConv, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_F32); + psArrayAdd(job->args, 1, footprintImage); 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 @@ -168,9 +188,9 @@ } else { psScalar *scalar = NULL; - scalar = job->args->data[5]; + scalar = job->args->data[8]; Next += scalar->data.S32; - scalar = job->args->data[6]; + scalar = job->args->data[9]; Npetro += scalar->data.S32; - scalar = job->args->data[7]; + scalar = job->args->data[10]; Nannuli += scalar->data.S32; } @@ -180,4 +200,5 @@ psFree (cellGroups); psFree(AnalysisRegion); + psFree(footprintImage); psLogMsg ("psphot", PS_LOG_WARN, "galaxy parameter analysis: %f sec for %d objects\n", psTimerMark ("psphot.galaxy"), Next); @@ -211,4 +232,5 @@ int Npetro = 0; int Nannuli = 0; + int Ngood = 0; // arguments: readout, sources, models, region, psfSize, maskVal, markVal @@ -218,5 +240,10 @@ psMetadata *recipe = job->args->data[3]; +#ifndef USE_FOOTPRINTS float skynoise = PS_SCALAR_VALUE(job->args->data[4],F32); +#endif + float fitNsigmaConv = PS_SCALAR_VALUE(job->args->data[5],F32); + float psfSize = PS_SCALAR_VALUE(job->args->data[6],F32); + psImage *footprintImage = job->args->data[7]; bool status; @@ -227,8 +254,18 @@ assert (maskVal); + // don't try and use bad models + pmModelStatus badModel = PM_MODEL_STATUS_NONE; + badModel |= PM_MODEL_STATUS_BADARGS; + badModel |= PM_MODEL_STATUS_OFFIMAGE; + badModel |= PM_MODEL_STATUS_NAN_CHISQ; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID; + badModel |= PM_MODEL_PCM_FAIL_GUESS; + pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC"); psVector *bgPixels = NULL; psVector *bgPixelsFlip = NULL; +#define dontMAKE_MARK_IMAGE 1 #ifdef MAKE_MARK_IMAGE psImage *markImage = NULL; @@ -239,4 +276,11 @@ #define MARK_BG 2 psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + +#define dontDUMP_IMAGES 1 +#ifdef DUMP_IMAGES + char *filebase = "/data/ippc18.0/home/bills/debug/psphot/images"; + char filename[80]; + psMetadata *dummyHeader = psMetadataAlloc(); +#endif // choose the sources of interest @@ -260,4 +304,7 @@ if (!source->extpars) SKIP (Nskip6); if (!source->modelFits) SKIP (Nskip7); + + // Flux cut to reduce sources to those that match SDSS + // if (source->moments->KronFlux < 2e5) continue; pmModel *sersicModel = NULL; @@ -270,4 +317,5 @@ } if (!sersicModel) SKIP (Nskip8); + if (sersicModel->flags & badModel) SKIP (Nskip8); Next ++; @@ -283,4 +331,5 @@ if (!isExtended || extModel != sersicModel) { + // XXX: I'm not doing this right. For now just look at objects for which source->modelEXT = psMemIncrRefCounter (sersicModel); source->type = PM_SOURCE_TYPE_EXTENDED; @@ -288,14 +337,10 @@ // 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 + if (source->modelEXT->isPCM) { + pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv); + } else { + pmSourceCacheModel (source, maskVal); + } + } psF32 *PAR = sersicModel->params->data.F32; @@ -306,6 +351,23 @@ psEllipseAxes axes = pmPSF_ModelToAxes (PAR, sersicModel->class->useReff); source->extpars->ghalfLightRadius = axes.major; + // and we will integrate out to 2 * half light radius psF32 R = 2 * axes.major; psF32 R2 = PS_SQR(R); + + 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); + // Make image at least 10 * half light radius + radius = PS_MAX (radius, 5*axes.major); + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 2*radius); + + +#ifdef DUMP_IMAGES + sprintf(filename, "%s/%04d-pix.fits", filebase, source->id); + + psphotSaveImage(dummyHeader, source->pixels, filename); +#endif + // Create the markImage and accumlate the Flux @@ -336,10 +398,98 @@ #endif +#ifdef USE_FOOTPRINTS + psS32 **vFootprint = footprintImage->data.S32; +#else // XXX use recipe value instead of hard coded 1.5 - psF32 threshold = skynoise * 1.5; + psF32 threshold = skynoise * 3; +#endif psF32 sumI = 0; + psS32 row0 = source->pixels->row0; + psS32 col0 = source->pixels->col0; 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; + psS32 yImage = row + row0; + 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 + + psS32 xImage = col + col0; + + psF32 xDiff = col - xCM; + + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 <= R2) { + sumI += pixel; + } + +#ifdef USE_FOOTPRINTS + if (vFootprint[yImage][xImage]) { +#else + if (fabs(pixel) >= threshold) { +#endif + 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; + +#ifdef USE_FOOTPRINTS + if (vFootprint[yFlip + row0][xFlip + col0] == 0) { +#else + if (fabs(mirrorPixel) < threshold) { +#endif + // This is a good background pixel with a good mirror Save the values + SET_MARK(MARK_BG); + psVectorAppend(bgPixels, pixel); + psVectorAppend(bgPixelsFlip, mirrorPixel); + } + } + } +#ifdef DUMP_IMAGES + sprintf(filename, "%s/%04d-mrk.fits", filebase, source->id); + + psphotSaveImage(dummyHeader, markImage, filename); +#endif + + // 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 + fprintf(stderr, "failed to subtract sersic model for source %d\n", source->id); + goto restoreModel; + } +#ifdef DUMP_IMAGES + sprintf(filename, "%s/%04d-res.fits", filebase, source->id); + + psphotSaveImage(dummyHeader, source->pixels, filename); +#endif + + // don't know if these pointers could change in the background subtraction so grab them again + 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; + if (fabs(yDiff) > R) continue; // y coordinate of mirror pixel int yFlip = yCM - yDiff; @@ -351,20 +501,8 @@ 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. + if (r2 > R2) continue; // x coordinate of mirror pixel @@ -376,51 +514,4 @@ 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); @@ -429,4 +520,11 @@ } } + + if (bgPixels->n < numPixels) { + fprintf(stderr, "could not find enough background pixels for source %d found %ld need %d\n", source->id, bgPixels->n, numPixels); + goto restoreModel; + } + Ngood++; + // and calculate Bt and Ba psF32 sumBt = 0; @@ -439,5 +537,9 @@ double rand = psRandomUniform(rng); long index = floor(rand * length); - if (index >= length) continue; + if (index >= length) { + fprintf(stderr, "random index %ld is larger than vector length %ld\n", + index, length); + continue; + } sumBt += fabs(vB[index] + vBFlip[index]); sumBa += fabs(vB[index] - vBFlip[index]); @@ -449,4 +551,11 @@ source->extpars->gS2 = source->extpars->gRT + source->extpars->gRA; +#ifdef PRINT_STUFF + if (Ngood % 40 == 1) { + printf(" id sumI rHL rad sumI gRT gRA sumRt sumBt numPixels numBGPixels\n"); + } + printf("%5d %5.2e %3.1f %10.1f %6.3f %6.3f %5.2e %5.2e %10d %6ld\n", + source->id, sumI, source->extpars->ghalfLightRadius, radius, source->extpars->gRT, source->extpars->gRA, sumRt, sumBt, numPixels, bgPixels->n); +#endif // Now put things back the way that they were. @@ -468,9 +577,13 @@ pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); } - } #ifdef MAKE_MARK_IMAGE psFree(markImage); +#endif +#ifdef DUMP_IMAGES + psFree(dummyHeader); + // psphotSaveImage raises an error + psErrorClear(); #endif psFree(bgPixels) @@ -481,11 +594,11 @@ // change the value of a scalar on the array (wrap this and put it in psArray.h) - scalar = job->args->data[5]; + scalar = job->args->data[8]; scalar->data.S32 = Next; - scalar = job->args->data[6]; + scalar = job->args->data[9]; scalar->data.S32 = Npetro; - scalar = job->args->data[7]; + scalar = job->args->data[10]; scalar->data.S32 = Nannuli;