Index: /branches/eam_branches/ipp-20130419/psModules/src/camera/pmReadoutFake.c =================================================================== --- /branches/eam_branches/ipp-20130419/psModules/src/camera/pmReadoutFake.c (revision 35550) +++ /branches/eam_branches/ipp-20130419/psModules/src/camera/pmReadoutFake.c (revision 35551) @@ -152,5 +152,5 @@ return false; } - + psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n", fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS], Index: /branches/eam_branches/ipp-20130419/psModules/src/detrend/pmDetrendDB.c =================================================================== --- /branches/eam_branches/ipp-20130419/psModules/src/detrend/pmDetrendDB.c (revision 35550) +++ /branches/eam_branches/ipp-20130419/psModules/src/detrend/pmDetrendDB.c (revision 35551) @@ -110,4 +110,5 @@ DETREND_STRING_CASE(VIDEODARK); DETREND_STRING_CASE(LINEARITY); + DETREND_STRING_CASE(AUXMASK); default: return NULL; Index: /branches/eam_branches/ipp-20130419/psModules/src/detrend/pmDetrendDB.h =================================================================== --- /branches/eam_branches/ipp-20130419/psModules/src/detrend/pmDetrendDB.h (revision 35550) +++ /branches/eam_branches/ipp-20130419/psModules/src/detrend/pmDetrendDB.h (revision 35551) @@ -39,4 +39,5 @@ PM_DETREND_TYPE_VIDEODARK, PM_DETREND_TYPE_LINEARITY, + PM_DETREND_TYPE_AUXMASK, } pmDetrendType; Index: /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmPSFEnvelope.c (revision 35550) +++ /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmPSFEnvelope.c (revision 35551) @@ -37,5 +37,5 @@ -//#define TESTING // Enable test output +// #define TESTING // Enable test output // #define PEAK_NORM // Normalise peaks? #define PEAK_FLUX 1.0e4 // Peak flux for each source @@ -190,6 +190,8 @@ psf->residuals = NULL; pmModelClassSetLimits(PM_MODEL_LIMITS_MODERATE); + psLogMsg("psModules",PS_LOG_INFO,"Matching Input %d",i); +#define CIRCULARIZE true if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, 0, xOffset, yOffset, psf, - NAN, radius, true, false)) { + NAN, radius, CIRCULARIZE, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout."); psFree(envelope); Index: /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmStack.c =================================================================== --- /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmStack.c (revision 35550) +++ /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmStack.c (revision 35551) @@ -41,6 +41,8 @@ /* #define TEST_X 5745 // x coordinate to examine */ /* #define TEST_Y 5331 // y coordinate to examine */ -#define TEST_X 25 -#define TEST_Y 25 +// #define TEST_X 972 +// #define TEST_Y 3213 +#define TEST_X 3289 +#define TEST_Y 4810 #define TEST_RADIUS 0.5 // Radius to examine # endif @@ -752,8 +754,12 @@ if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { for (int i = 0; i < numGood; i++) { - fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n", - i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i], - addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i], - pixelSuspects->data.U8[i], badMaskBits, suspectMaskBits, *badMask, *goodMask); + fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%g) %g %f %d %x %x -> %x %x\n", + i, x, y, pixelSources->data.U16[i], + pixelData->data.F32[i], pixelVariances->data.F32[i], + addVariance->data.F32[i], + pixelWeights->data.F32[i], pixelExps->data.F32[i], + pixelSuspects->data.U8[i], + badMaskBits, suspectMaskBits, + *badMask, *goodMask); } } Index: /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmSubtraction.c (revision 35550) +++ /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmSubtraction.c (revision 35551) @@ -37,4 +37,5 @@ #define USE_KERNEL_ERR // Use kernel error image? #define NUM_COVAR_POS 5 // Number of positions for covariance calculation +#define USE_LOGFIT_REJECT // XXX we need to pass these fwhm values elsewhere. These should go on one of the structure, but @@ -966,4 +967,78 @@ psPolynomial1D *model = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2); +#ifdef USE_LOGFIT_REJECT + // CZW: Since flux and chisq can range over many orders of magnitude, use a log-log fit to calculate the + // RMS scatter. This should be more robust against outliers that can impact the fit quality. + psVector *logchisq = psVectorAlloc(match->chisq->n,PS_TYPE_F32); + psVector *logflux = psVectorAlloc(match->fluxes->n,PS_TYPE_F32); + for (int qq = 0; qq < match->chisq->n; qq++) { + if ((match->chisq->data.F32[qq] > 0)&& + (match->fluxes->data.F32[qq] > 0)) { + logchisq->data.F32[qq] = log10(match->chisq->data.F32[qq]); + logflux->data.F32[qq] = log10(match->fluxes->data.F32[qq]); + } + else { + // Ignore negative values. Maybe we could do an offset or something, but not today. + match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[qq] |= PM_SUBTRACTION_STAMP_REJECTED; + } + // if (!(match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[qq] & 0xff)) { + // psLogMsg("psModules.imcombine", PS_LOG_INFO, "RRRRRR: %d %g %g %g %g\n", + // qq,match->chisq->data.F32[qq],match->fluxes->data.F32[qq], + // 0.0,0.0); + // } + + } + bool result = psVectorClipFitPolynomial1D(model, stats, match->stampMask, 0xff, logchisq, NULL, logflux); + if (!result) { + psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations."); + psFree(model); + psFree(stats); + return -1; + } + + // However, since we've done a log-log fit, we can't rely on the statistics in the stats object + // To accurately represent the distribution. I've tried to massage them (as dlog10(X) = abs(1/(X log(10))) dX), + // but ended up deciding to just manually calculate the residual, and then pass that to a robust stats object. + psStats *residStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psVector *residData = psVectorAlloc(match->chisq->n,PS_TYPE_F32); + int counter = 0; + for (int qq = 0; qq < match->chisq->n; qq++) { + if (!(match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[qq] & 0xff)) { + double vv = pow(10, + model->coeff[0] + + model->coeff[1] * log10(match->fluxes->data.F32[qq]) + + model->coeff[2] * pow(log10(match->fluxes->data.F32[qq]),2)); + residData->data.F32[qq] = match->chisq->data.F32[qq] - vv; + counter++; + // psLogMsg("psModules.imcombine", PS_LOG_INFO, "SSSSS: %d %g %g %g %g\n", + // qq,match->chisq->data.F32[qq],match->fluxes->data.F32[qq], + // vv,residData->data.F32[qq]); + } + } + psVectorStats(residStats,residData,NULL,match->stampMask,0xff); + if (isnan(residStats->robustMedian)) { + psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations."); + psFree(model); + psFree(stats); + psFree(logchisq); + psFree(logflux); + psFree(residStats); + psFree(residData); + + return -1; + } + + kernels->mean = residStats->robustMedian; + kernels->rms = residStats->robustStdev; + kernels->numStamps = counter; + + psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %f +/- %f\n", residStats->robustMedian,residStats->robustStdev); + psFree(logchisq); + psFree(logflux); + psFree(residStats); + psFree(residData); + +#else + // CZW: Otherwise, use the original linear fit code. bool result = psVectorClipFitPolynomial1D(model, stats, match->stampMask, 0xff, match->chisq, NULL, match->fluxes); if (!result) { @@ -979,11 +1054,11 @@ return -1; } - kernels->mean = stats->sampleMean; - kernels->rms = stats->sampleStdev; + kernels->rms = stats->sampleStdev; kernels->numStamps = stats->clippedNvalues; + psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %f +/- %f\n", stats->sampleMean, stats->sampleStdev); +#endif psLogMsg ("pmPSFtry", 4, "chisq vs flux model: %e + %e flux + %e flux^2\n", model->coeff[0], model->coeff[1], model->coeff[2]); - psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %f +/- %f\n", stats->sampleMean, stats->sampleStdev); psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf", kernels->numStamps, kernels->mean, kernels->rms); Index: /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmSubtractionEquation.c (revision 35550) +++ /branches/eam_branches/ipp-20130419/psModules/src/imcombine/pmSubtractionEquation.c (revision 35551) @@ -751,5 +751,5 @@ stamps->normValue2 = stats->robustMedian; - psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2); + psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f) %ld\n", stamps->normValue, stamps->normValue2,norms->n); psFree(stats); @@ -790,5 +790,7 @@ } } - + // psLogMsg("psModules.imcombine",PS_LOG_INFO, "stampNorm: %f %f %d %d %d %f %f", + // stamp->x,stamp->y,footprint,normWindow1,normWindow2, + // normI1,normI2); stamp->norm = normI2 / normI1; stamp->normI1 = normI1; Index: /branches/eam_branches/ipp-20130419/psModules/src/objects/pmPSF.c =================================================================== --- /branches/eam_branches/ipp-20130419/psModules/src/objects/pmPSF.c (revision 35550) +++ /branches/eam_branches/ipp-20130419/psModules/src/objects/pmPSF.c (revision 35551) @@ -342,5 +342,5 @@ PS_ASSERT_PTR_NON_NULL(modelPar, axes); - bool useReff = true; + bool useReff = false; useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC")); useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV")); @@ -386,13 +386,13 @@ psEllipseShape shape = psEllipseAxesToShape (axes); - bool useReff = true; - useReff |= pmModelClassGetType ("PS_MODEL_SERSIC"); - useReff |= pmModelClassGetType ("PS_MODEL_DEV"); - useReff |= pmModelClassGetType ("PS_MODEL_EXP"); + bool useReff = false; + useReff |= ( type == pmModelClassGetType ("PS_MODEL_SERSIC")); + useReff |= ( type == pmModelClassGetType ("PS_MODEL_DEV")); + useReff |= ( type == pmModelClassGetType ("PS_MODEL_EXP")); if (useReff) { modelPar[PM_PAR_SXX] = shape.sx; modelPar[PM_PAR_SYY] = shape.sy; - modelPar[PM_PAR_SXY] = shape.sxy; // XXX NEED factor of 2 here for correct angle conversion + modelPar[PM_PAR_SXY] = shape.sxy * 2.0; // XXX NEED factor of 2 here for correct angle conversion } else { modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2;