Index: /branches/pap/psphot/src/psphotFake.c =================================================================== --- /branches/pap/psphot/src/psphotFake.c (revision 25311) +++ /branches/pap/psphot/src/psphotFake.c (revision 25312) @@ -4,5 +4,5 @@ PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models -#define TESTING +//#define TESTING @@ -12,5 +12,7 @@ // non-Gaussian PSF) to the limiting flux in the un-smoothed original image. static bool fakeLimit(float *magLim, // Limiting magntiude, to return - int *radius, // Radius for fake sources + int *radius, // Radius for fake sources, to return + float *minFlux, // Minimum flux for fake sources, to return + float *norm, // Normalisation of PSF (conversion: peak --> integrated flux) const pmReadout *ro, // Readout of interest const pmPSF *psf, // Point-spread function @@ -25,4 +27,6 @@ assert(magLim); assert(radius); + assert(minFlux); + assert(norm); psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing @@ -53,5 +57,5 @@ return false; } - float norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0 + *norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0 psFree(normModel); @@ -65,5 +69,5 @@ // I_smooth = Flux / 2*norm. float peakLim = thresh * sqrtf(meanVar * factor); // Limiting peak value in smoothed image - float fluxLim = 2.0 * norm * peakLim; // Limiting flux in original + float fluxLim = 2.0 * *norm * peakLim; // Limiting flux in original *magLim = -2.5 * log10f(fluxLim); psTrace("psphot.fake", 1, "Limiting peak: %f\n", peakLim); @@ -72,4 +76,6 @@ *radius = smoothSigma * smoothNsigma; + + *minFlux = 0.1 * sqrtf(meanVar); return true; @@ -83,5 +89,6 @@ int numSources, // Number of fake sources for each bin float refMag, // Reference magnitude - int radius // Radius for fake sources + int radius, // Radius for fake sources + float minFlux // Minimum flux for fake sources ) { @@ -122,5 +129,5 @@ pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout if (!pmReadoutFakeFromVectors(fakeRO, numCols, numRows, xAll, yAll, magAll, - NULL, NULL, psf, NAN, radius, false, true)) { + NULL, NULL, psf, minFlux, radius, false, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image"); psFree(fakeRO); @@ -191,5 +198,5 @@ - +#if defined(TESTING) && 0 { psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); @@ -205,5 +212,5 @@ psFree(rng); } - +#endif @@ -211,5 +218,8 @@ float magLim; // Guess at limiting magnitude int radius; // Radius for fake sources - if (!fakeLimit(&magLim, &radius, readout, psf, thresh, smoothSigma, smoothNsigma, maskVal)) { + float minFlux; // Minimum flux for fake sources + float norm; // Normalisation of PSF + if (!fakeLimit(&magLim, &radius, &minFlux, &norm, readout, + psf, thresh, smoothSigma, smoothNsigma, maskVal)) { psError(PS_ERR_UNKNOWN, false, "Unable to determine limits for image"); return false; @@ -223,5 +233,5 @@ psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row - if (!fakeGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius)) { + if (!fakeGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius, minFlux)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources"); psFree(xFake); @@ -246,5 +256,5 @@ #ifdef TESTING - psphotSaveImage(NULL, readout->mask, "fake_sig.fits"); + psphotSaveImage(NULL, significance, "fake_sig.fits"); #endif @@ -253,77 +263,99 @@ int numBins = magOffsets->n; // Number of bins psVector *count = psVectorAlloc(numBins, PS_TYPE_S32); // Number of sources found in each bin - psArray *fakeSources = psArrayAllocEmpty(numSources * numBins); // Fake sources - psVector *fakeEnd = psVectorAllocEmpty(numSources * numBins, PS_TYPE_S32); // End index of each bin - long numFound = 0; // Number of sources found + psArray *fakeSources = psArrayAlloc(numSources); // Fake sources in each bin + psArray *fakeSourcesAll = psArrayAllocEmpty(numSources * numBins); // All fake sources for (int i = 0; i < numBins; i++) { - int num = 0; // Number found - + int numFound = 0; // Number found + + // Determine extraction size + float mag = magLim + magOffsets->data.F32[i]; // Magnitude for bin + float peak = powf(10.0, -0.4 * mag) / norm; // Peak flux + pmModel *model = pmModelFromPSFforXY(psf, numCols / 2.0, numRows / 2.0, peak); // Model for source + if (!model || (model->flags & MODEL_MASK)) { + psError(PS_ERR_UNKNOWN, true, "Unable to generate model for bin %d", i); + psFree(model); + return false; + } + float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source + psFree(model); + + psArray *sources = psArrayAllocEmpty(numSources); // Sources in this bin for (int j = 0; j < numSources; j++) { // Coordinates of interest - float x = PS_MIN(xFake->data.F32[i][j] + 0.5, numCols - 1); - float y = PS_MIN(yFake->data.F32[i][j] + 0.5, numRows - 1); - int xPix = x, yPix = y; // Pixel coordinates + float x = xFake->data.F32[i][j]; + float y = yFake->data.F32[i][j]; + int xPix = PS_MIN(x + 0.5, numCols - 1), yPix = PS_MIN(y + 0.5, numRows - 1); // Pixel coordinates float sig = significance->data.F32[yPix][xPix]; // Significance of pixel if (sig > thresh) { pmSource *source = pmSourceAlloc(); // Fake source - source->peak = pmPeakAlloc(x, y, NAN, PM_PEAK_LONE); - if (!pmSourceDefinePixels(source, readout, x, y, radius)) { + source->peak = pmPeakAlloc(x, y, sig, PM_PEAK_LONE); + if (!pmSourceDefinePixels(source, readout, x, y, sourceRadius)) { psErrorClear(); continue; } + source->peak->xf = x; + source->peak->yf = y; source->modelPSF = pmModelFromPSFforXY(psf, x, y, 2.0 * sqrtf(sig)); source->type = PM_SOURCE_TYPE_STAR; - num++; - - fakeSources->data[numFound] = source; numFound++; + psArrayAdd(sources, sources->n, source); + psArrayAdd(fakeSourcesAll, fakeSourcesAll->n, source); + psFree(source); } } - count->data.S32[i] = num; - fakeEnd->data.S32[i] = numFound; - } - fakeSources->n = numFound; - + fakeSources->data[i] = sources; + count->data.S32[i] = numFound; + } psFree(xFake); psFree(yFake); psFree(significance); - if (!psphotFitSourcesLinear(readout, fakeSources, recipe, psf, true)) { + if (!psphotFitSourcesLinear(readout, fakeSourcesAll, recipe, psf, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); psFree(fakeSources); - psFree(fakeEnd); psFree(count); return false; } - if (!psphotMagnitudes(config, readout, view, fakeSources, psf)) { + + // Disable aperture corrections; casting away const! + pmTrend2D *apTrend = psf->ApTrend; // Aperture trend + ((pmPSF*)psf)->ApTrend = NULL; + + if (!psphotMagnitudes(config, readout, view, fakeSourcesAll, psf)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources."); psFree(fakeSources); - psFree(fakeEnd); psFree(count); - return false; - } - + ((pmPSF*)psf)->ApTrend = apTrend; // Casting away const! + return false; + } + psFree(fakeSourcesAll); + + // Re-enable aperture corrections; casting away const! + ((pmPSF*)psf)->ApTrend = apTrend; + + psVector *magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean difference in magnitude for each bin + psVector *magDiffStdev = psVectorAlloc(numBins, PS_TYPE_F32); // Stdev of diff in magnitude for each bin psVector *magErrMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean error in magnitude for each bin - psVector *magErrStdev = psVectorAlloc(numBins, PS_TYPE_F32); // Stdev of error in magnitude for each bin - psVector *magErrErrMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean error in magnitude for each bin + psVector *magDiff = psVectorAlloc(numSources, PS_TYPE_F32); // Magnitude differences psVector *magErr = psVectorAlloc(numSources, PS_TYPE_F32); // Magnitude errors - psVector *magErrErr = psVectorAlloc(numSources, PS_TYPE_F32); // Error in magnitude error psVector *magMask = psVectorAlloc(numSources, PS_TYPE_VECTOR_MASK); // Mask for magnitude errors - psStats *magStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics - psVector *robustMedian = psVectorAlloc(numBins, PS_TYPE_F32); - psVector *robustStdev = psVectorAlloc(numBins, PS_TYPE_F32); - for (int i = 0, index = 0; i < numBins; i++) { + psStats *magStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | + PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics + for (int i = 0; i < numBins; i++) { psStatsInit(magStats); psVectorInit(magMask, 0); +#ifdef TESTING psString name = NULL; psStringAppend(&name, "fake_%d.dat", i); FILE *file = fopen(name, "w"); + psFree(name); +#endif float magRef = magLim + magOffsets->data.F32[i]; // Reference magnitude - int j; // Index - for (j = 0; index < fakeEnd->data.S32[i]; j++, index++) { - pmSource *source = fakeSources->data[index]; // Source of interest + psArray *sources = fakeSources->data[i]; // Sources in bin + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; // Source of interest if (!source || !isfinite(source->psfMag)) { magMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF; @@ -331,34 +363,51 @@ } - fprintf(file, "%f %f %f %f %f\n", source->peak->xf, source->peak->yf, +#ifdef TESTING + fprintf(file, "%f %f %f %f %f %f %f\n", source->peak->xf, source->peak->yf, + source->modelPSF->params->data.F32[PM_PAR_XPOS], + source->modelPSF->params->data.F32[PM_PAR_YPOS], magRef, source->psfMag, source->errMag); - magErr->data.F32[j] = source->psfMag - magRef; - magErrErr->data.F32[j] = source->errMag; - } - magErr->n = j; - magMask->n = j; +#endif + magDiff->data.F32[j] = source->psfMag - magRef; + magErr->data.F32[j] = source->errMag; + } + magDiff->n = sources->n; + magMask->n = sources->n; + magErr->n = sources->n; + if (!psVectorStats(magStats, magDiff, NULL, magMask, 0xFF)) { + // Probably because we don't have enough sources + psErrorClear(); + } + + if (isfinite(magStats->robustMedian) && isfinite(magStats->robustStdev)) { + magDiffMean->data.F32[i] = magStats->robustMedian; + magDiffStdev->data.F32[i] = magStats->robustStdev; + } else { + magDiffMean->data.F32[i] = magStats->sampleMean; + magDiffStdev->data.F32[i] = magStats->sampleStdev; + } + if (!psVectorStats(magStats, magErr, NULL, magMask, 0xFF)) { // Probably because we don't have enough sources psErrorClear(); } - magErrMean->data.F32[i] = magStats->sampleMean; - magErrStdev->data.F32[i] = magStats->sampleStdev; - robustMedian->data.F32[i] = magStats->robustMedian; - robustStdev->data.F32[i] = magStats->robustStdev; - if (!psVectorStats(magStats, magErrErr, NULL, magMask, 0xFF)) { - // Probably because we don't have enough sources - psErrorClear(); - } - magErrErrMean->data.F32[i] = magStats->sampleMean; + + if (isfinite(magStats->robustMedian)) { + magErrMean->data.F32[i] = magStats->robustMedian; + } else { + magErrMean->data.F32[i] = magStats->sampleMean; + } + +#ifdef TESTING fclose(file); - } - +#endif + } psFree(magStats); + psFree(magDiff); psFree(magMask); psFree(magErr); psFree(fakeSources); - psFree(fakeEnd); // XXX How do we get the results out? @@ -372,23 +421,17 @@ psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.COUNTS", PS_META_REPLACE, "Number of sources retrieved", count); - psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.ERR", PS_META_REPLACE, - "Mean magnitude differences", magErrMean); - psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.RMS", PS_META_REPLACE, - "Stdev in magnitude differences", magErrStdev); - psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.RM", PS_META_REPLACE, - "Robust median magnitude differences", robustMedian); - psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.RS", PS_META_REPLACE, - "Robust stdev in magnitude differences", robustStdev); - psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.MEANERR", PS_META_REPLACE, - "Mean error in magnitude differences", magErrErrMean); + psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.DIFF.MEAN", PS_META_REPLACE, + "Mean magnitude differences", magDiffMean); + psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.DIFF.STDEV", PS_META_REPLACE, + "Stdev of magnitude differences", magDiffStdev); + psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.ERR.MEAN", PS_META_REPLACE, + "Mean error in magnitude differences", magErrMean); psMetadataConfigWrite(stats, "fake.stats"); psFree(stats); psFree(count); + psFree(magDiffMean); + psFree(magDiffStdev); psFree(magErrMean); - psFree(magErrStdev); - - psFree(robustMedian); - psFree(robustStdev); return true;