Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 5083) +++ /trunk/psphot/Makefile (revision 5084) @@ -15,5 +15,6 @@ # LIBS = -lpsmodule $(LPSLIB) LIBS = $(LPSLIB) -CFLAGS = $(INCS) -std=c99 -Wall -Werror -g +# CFLAGS = $(INCS) -std=c99 -Wall -Werror -g +CFLAGS = $(INCS) -std=c99 -g LFLAGS = $(LIBS) Index: /trunk/psphot/src/pmPSFtry.c =================================================================== --- /trunk/psphot/src/pmPSFtry.c (revision 5083) +++ /trunk/psphot/src/pmPSFtry.c (revision 5084) @@ -161,5 +161,5 @@ // XXX this function wants aperture radius for pmSourcePhotometry - pmPSFtryMetric (try, RADIUS); + pmPSFtryMetric_Alt (try, RADIUS); psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n", modelName, try->psf->ApResid, try->psf->dApResid, try->psf->skyBias); @@ -167,5 +167,4 @@ return (try); } - bool pmPSFtryMetric (pmPSFtry *try, float RADIUS) { @@ -195,4 +194,6 @@ } fclose (f); + + // XXX EAM : try 3hi/1lo sigma clipping on the rflux v dap fit // find min and max of (1/flux): @@ -297,2 +298,69 @@ return true; } + +bool pmPSFtryMetric_Alt (pmPSFtry *try, float RADIUS) { + + // the measured (aperture - fit) magnitudes (dA == try->metric) + // depend on both the true ap-fit (dAo) and the bias in the sky measurement: + // dA = dAo + dsky/flux + // where flux is the flux of the star + // we fit this trend to find the infinite flux aperture correction (dAo), + // the nominal sky bias (dsky), and the error on dAo + // the values of dA are contaminated by stars with close neighbors in the aperture + // we use an outlier rejection to avoid this bias + + // rflux = ten(0.4*fitMag); + psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64); + for (int i = 0; i < try->sources->n; i++) { + if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; + rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]); + } + + // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + + // XXX EAM + stats->min = 1.0; + stats->max = 3.0; + stats->clipIter = 3; + + // linear fit to rfBin, daBin + psPolynomial1D *poly = Polynomial1DAlloc_EAM (PS_POLYNOMIAL_ORD, 1); + + // XXX EAM : this is the intended API (cycle 7? cycle 8?) + poly = VectorClipFitPolynomial1D_EAM (poly, stats, try->mask, PSFTRY_MASK_ALL, try->metric, NULL, rflux); + + fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); + + // these stats should come back from the fit above... + // psVector *daFit = Polynomial1DEvalVector_EAM (poly, rflux); + // psVector *daResid = (psVector *) psBinaryOp (NULL, (void *) try->metric, "-", (void *) daFit); + + // stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV); + // stats->clipIter = 3; + // stats->clipSigma = 3; + + // stats = psVectorStats (stats, daResid, NULL, maskB, 1); + + try->psf->ApResid = poly->coeff[0]; + try->psf->dApResid = stats->sampleStdev; + try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); + + FILE *f; + f = fopen ("apresid.dat", "w"); + if (f == NULL) psAbort ("pmPSFtry", "can't open output file"); + + for (int i = 0; i < try->sources->n; i++) { + fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i], try->mask->data.U8[i]); + } + fclose (f); + + psFree (rflux); + psFree (poly); + psFree (stats); + + // psFree (daFit); + // psFree (daResid); + + return true; +} Index: /trunk/psphot/src/pmPSFtry.h =================================================================== --- /trunk/psphot/src/pmPSFtry.h (revision 5083) +++ /trunk/psphot/src/pmPSFtry.h (revision 5084) @@ -27,4 +27,5 @@ pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float radius); bool pmPSFtryMetric (pmPSFtry *try, float RADIUS); +bool pmPSFtryMetric_Alt (pmPSFtry *try, float RADIUS); # endif Index: /trunk/psphot/src/psPolynomials.c =================================================================== --- /trunk/psphot/src/psPolynomials.c (revision 5083) +++ /trunk/psphot/src/psPolynomials.c (revision 5084) @@ -507,9 +507,13 @@ psVector *zFit = NULL; psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64); - - // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats: + float minClipSigma = stats->min; + float maxClipSigma = stats->max; + float minClipValue; + float maxClipValue; + + // XXX EAM : use SAMPLE_MEDIAN and SAMPLE_STDEV for stats: stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - for (int N = 0; N < 3; N++) { + for (int N = 0; N < stats->clipIter; N++) { int Nkeep = 0; @@ -519,5 +523,8 @@ stats = psVectorStats (stats, zResid, NULL, mask, maskValue); - psTrace (".psphot.RobustFit", 4, "residual stats for robust fit: %g +/- %g\n", + minClipValue = -1*fabs(minClipSigma)*stats->sampleStdev; + maxClipValue = +1*fabs(maxClipSigma)*stats->sampleStdev; + + psTrace (".psphot.VectorClipFit", 4, "residual stats for robust fit: %g +/- %g\n", stats->sampleMedian, stats->sampleStdev); @@ -527,5 +534,9 @@ for (int i = 0; i < zResid->n; i++) { if (mask->data.U8[i]) continue; - if (fabs(zResid->data.F64[i] - stats->sampleMedian) > 2*stats->sampleStdev) { + if (zResid->data.F64[i] - stats->sampleMedian > maxClipValue) { + mask->data.U8[i] |= 0x01; + continue; + } + if (zResid->data.F64[i] - stats->sampleMedian < minClipValue) { mask->data.U8[i] |= 0x01; continue; @@ -534,5 +545,5 @@ } - psTrace (".psphot.RobustFit", 4, "keeping %d of %d pts for fit\n", + psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n", Nkeep, x->n); Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 5083) +++ /trunk/psphot/src/psphot.h (revision 5084) @@ -39,9 +39,9 @@ // output functions bool pmSourcesWriteText (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -int pmSourcesDophotType (pmSource *source); +bool pmSourcesWriteOBJ (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteCMP (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteCMF (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); + bool pmPeaksWriteText (psArray *sources, char *filename); bool pmMomentsWriteText (psArray *sources, char *filename); @@ -49,2 +49,4 @@ bool pmModelWriteFLTs (psArray *sources, char *filename); bool pmModelWriteNULLs (psArray *sources, char *filename); + +int pmSourcesDophotType (pmSource *source); Index: /trunk/psphot/src/psphotArguments.c =================================================================== --- /trunk/psphot/src/psphotArguments.c (revision 5083) +++ /trunk/psphot/src/psphotArguments.c (revision 5084) @@ -38,4 +38,12 @@ } + // optional output residual image - add to config + char *photcode = NULL; + if ((N = psArgumentGet (*argc, argv, "-photcode"))) { + psArgumentRemove (N, argc, argv); + photcode = psStringCopy (argv[N]); + psArgumentRemove (N, argc, argv); + } + if (*argc != 4) usage (); @@ -60,4 +68,7 @@ psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid); } + if (photcode != NULL) { + psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", photcode); + } return (config); } @@ -70,4 +81,5 @@ fprintf (stderr, " -weight (filename)\n"); fprintf (stderr, " -resid (filename)\n"); + fprintf (stderr, " -photcode (photcode)\n"); exit (2); } Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 5083) +++ /trunk/psphot/src/psphotOutput.c (revision 5084) @@ -26,5 +26,5 @@ if (!strcasecmp (outputMode, "OBJ")) { - pmSourcesWriteOBJ (imdata, outputFile, sources, psf, sky); + pmSourcesWriteOBJ (imdata, config, outputFile, sources, psf, sky); return; } @@ -36,10 +36,10 @@ if (!strcasecmp (outputMode, "CMP")) { - pmSourcesWriteCMP (imdata, outputFile, sources, psf, sky); + pmSourcesWriteCMP (imdata, config, outputFile, sources, psf, sky); return; } if (!strcasecmp (outputMode, "CMF")) { - pmSourcesWriteCMF (imdata, outputFile, sources, psf, sky); + pmSourcesWriteCMF (imdata, config, outputFile, sources, psf, sky); return; } @@ -68,5 +68,5 @@ // dophot-style output list with fixed line width -bool pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { +bool pmSourcesWriteOBJ (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { int i, type, status; @@ -74,4 +74,6 @@ psF32 *PAR, *dPAR; float sky, dmag, apMag, fitMag; + float x, y, rflux; + bool result; psLine *line = psLineAlloc (104); // 104 is dophot-defined line length @@ -82,4 +84,6 @@ return false; } + + float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS"); // write sources with models first @@ -296,14 +300,24 @@ // elixir-style pseudo FITS table (header + ascii list) -bool pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { +bool pmSourcesWriteCMP (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { int i, type, status; pmModel *model; + psMetadataItem *mdi; psF32 *PAR, *dPAR; - float sky, dmag, apMag, fitMag, lsky, Theta; + float sky, dmag, apMag, fitMag, lsky; + float x, y, rflux; + bool result; // create file, write-out header unlink (filename); psFits *fits = psFitsAlloc (filename); + + // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0 + mdi = psMetadataLookup (imdata->header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_META_S32; + + // psMetadataAdd (imdata->header, PS_LIST_HEAD, "NAXIS", PS_META_S32 | PS_META_REPLACE, "", 0); psFitsWriteHeader (imdata->header, fits); psFree (fits); @@ -317,5 +331,16 @@ fseek (f, 0, SEEK_END); - psLine *line = psLineAlloc (66); // 66 is imclean-defined line length + psLine *line = psLineAlloc (67); // 66 is imclean-defined line length + + float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS"); + + // XXX EAM : need to add these concepts to the header (from psf and config) + // fits_modify (header, "ZERO_PT", "%lf", 1, ZERO_POINT); + // fits_modify (header, "FWHM_X", "%lf", 1, FWHMx); + // fits_modify (header, "FWHM_Y", "%lf", 1, FWHMy); + // fits_modify (header, "ANGLE", "%lf", 1, angle); + // fits_modify (header, "FSATUR", "%lf", 1, saturate); // can use 0.0 + // fits_modify (header, "FLIMIT", "%lf", 1, complete); // can use 0.0 + // fits_modify (header, "NSTARS", "%d", 1, N); // write sources with models first @@ -404,5 +429,5 @@ psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' - psLineAdd (line, "%5.1f\n", Theta); // should be theta + psLineAdd (line, "%5.1f\n", 0); // should be theta fwrite (line->line, 1, line->Nline, f); } @@ -414,15 +439,29 @@ // this format consists of a header derived from the image header // followed by a zero-size matrix, followed by the table data -bool pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { +bool pmSourcesWriteCMF (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { psArray *table; + psMetadataItem *mdi; psMetadata *row; psMetadata *theader; + int i, type, status; + pmModel *model; + psF32 *PAR, *dPAR; + float sky, dmag, apMag, fitMag; + float x, y, rflux, lsky; + bool result; // set NAXIS to 0 - psMetadataAdd (header, PS_LIST_HEAD, "NAXIS", PS_META_S32 | PS_META_REPLACE, "", 0); + fprintf (stderr, "setting naxis\n"); + mdi = psMetadataLookup (imdata->header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_META_S32; + + // psMetadataAdd (imdata->header, PS_LIST_HEAD, "NAXIS", PS_META_S32 | PS_META_REPLACE, "", 0); table = psArrayAlloc (sources->n); sources->n = 0; + + float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS"); for (i = 0; i < sources->n; i++) { @@ -499,16 +538,16 @@ lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); - table = psMetdataAlloc (); - psMetadataAdd (table, PS_LIST_HEAD, "X_PIX", PS_META_F32, "", PAR[2]); - psMetadataAdd (table, PS_LIST_HEAD, "Y_PIX", PS_META_F32, "", PAR[3]); - psMetadataAdd (table, PS_LIST_HEAD, "MAG_RAW", PS_META_F32, "", fitMag + 25.0); - psMetadataAdd (table, PS_LIST_HEAD, "MAG_ERR", PS_META_F32, "", (int)(1000*dmag)); - psMetadataAdd (table, PS_LIST_HEAD, "MAG_GAL", PS_META_F32, "", type); - psMetadataAdd (table, PS_LIST_HEAD, "MAG_AP", PS_META_F32, "", lsky); - psMetadataAdd (table, PS_LIST_HEAD, "LOG_SKY", PS_META_F32, "", 32.0); // should be 'Mgal' - psMetadataAdd (table, PS_LIST_HEAD, "FWHM_X", PS_META_F32, "", apMag + 25.0); // 25.0 should come from config - psMetadataAdd (table, PS_LIST_HEAD, "FWHM_Y", PS_META_F32, "", PAR[4]); // should be 'FHWM x' - psMetadataAdd (table, PS_LIST_HEAD, "THETA", PS_META_F32, "", PAR[5]); // should be 'FHWM y' - psMetadataAdd (table, PS_LIST_HEAD, "DOPHOT", PS_META_S8, "", Theta); // should be theta + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_HEAD, "X_PIX", PS_META_F32, "", PAR[2]); + psMetadataAdd (row, PS_LIST_HEAD, "Y_PIX", PS_META_F32, "", PAR[3]); + psMetadataAdd (row, PS_LIST_HEAD, "MAG_RAW", PS_META_F32, "", fitMag + 25.0); + psMetadataAdd (row, PS_LIST_HEAD, "MAG_ERR", PS_META_F32, "", (int)(1000*dmag)); + psMetadataAdd (row, PS_LIST_HEAD, "MAG_GAL", PS_META_F32, "", type); + psMetadataAdd (row, PS_LIST_HEAD, "MAG_AP", PS_META_F32, "", lsky); + psMetadataAdd (row, PS_LIST_HEAD, "LOG_SKY", PS_META_F32, "", 32.0); // should be 'Mgal' + psMetadataAdd (row, PS_LIST_HEAD, "FWHM_X", PS_META_F32, "", apMag + 25.0); // 25.0 should come from config + psMetadataAdd (row, PS_LIST_HEAD, "FWHM_Y", PS_META_F32, "", PAR[4]); // should be 'FHWM x' + psMetadataAdd (row, PS_LIST_HEAD, "THETA", PS_META_F32, "", PAR[5]); // should be 'FHWM y' + psMetadataAdd (row, PS_LIST_HEAD, "DOPHOT", PS_TYPE_S8, "", 0); // should be theta // psMetadataAdd (header, PS_LIST_HEAD, "DUMMY", PS_META_STR, "", NULL); @@ -517,10 +556,15 @@ } - image = psImageAlloc (0, 0, PS_TYPE_F32); - theader = psFitsHeaderFromTable (NULL, table, "SMPFILE"); + // psAbort ("psphotOutput", "FITS table output not finished"); + + psImage *image = psImageAlloc (1, 1, PS_TYPE_F32); + // XXX missing: theader = psFitsHeaderFromTable (NULL, table, "SMPFILE"); + + theader = psMetadataAlloc (); + psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_META_STR, "extension name", "SMPFILE"); unlink (filename); psFits *fits = psFitsAlloc (filename); - psFitsWriteImage (fits, header, image, 0, NULL); + psFitsWriteImage (fits, imdata->header, image, 0); psFitsWriteTable (fits, theader, table); psFree (fits);