Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 4640) +++ /trunk/psphot/Makefile (revision 4641) @@ -31,14 +31,10 @@ $(SRC)/apply_psf_model.$(ARCH).o \ $(SRC)/test_psf_scatter.$(ARCH).o \ -$(SRC)/mark_psf_sources.$(ARCH).o \ -$(SRC)/subtract_psf_sources.$(ARCH).o \ $(SRC)/mark_psf_source.$(ARCH).o \ $(SRC)/subtract_psf_source.$(ARCH).o \ $(SRC)/by_SN.$(ARCH).o \ $(SRC)/fit_galaxies.$(ARCH).o \ -$(SRC)/subtract_galaxies.$(ARCH).o \ $(SRC)/LocalSky.$(ARCH).o \ -$(SRC)/basic_classes.$(ARCH).o \ -$(SRC)/find_defects.$(ARCH).o +$(SRC)/basic_classes.$(ARCH).o FITSOURCE = \ Index: unk/psphot/src/find_defects.c =================================================================== --- /trunk/psphot/src/find_defects.c (revision 4640) +++ (revision ) @@ -1,127 +1,0 @@ -# include "psphot.h" -# define NSIGMA 3.0 - -bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf) -{ - - bool status; - float gradx, grady, gradlim; - float N0, N1, peak, sky; - int x, y, Nzap; - - psTimerStart ("psphot"); - - // measure the x and y max gradients - float RADIUS = psMetadataLookupF32 (&status, config, "FIT_RADIUS"); - int dN = 2*RADIUS + 1; - - psImage *flux = psImageAlloc (dN, dN, PS_TYPE_F32); - - // create an object located at the image center - psModel *model = psModelAlloc (psf->type); - model->params->data.F32[0] = 0; - model->params->data.F32[1] = 1; - // XXX : temporary hack : center of 2k x 4k image - model->params->data.F32[2] = 1024; - model->params->data.F32[3] = 2048; - - // fill in the PSF parameters - psModelFromPSFFunc modelFromPSFFunc = psModelFromPSFFunc_GetFunction (psf->type); - modelFromPSFFunc (model, model, psf); - - // generate the image flux - pmSourceAddModel (flux, NULL, model, true); - DumpImage (flux, "psf.fits"); - - fprintf (stderr, "sx: %f, sy: %f\n", model->params->data.F32[4], model->params->data.F32[5]); - - // measure max dfdx - gradx = 0; - for (int iy = 0; iy < flux->numRows; iy++) { - for (int ix = 0; ix < flux->numCols - 1; ix++) { - gradx = PS_MAX (gradx, fabs(flux->data.F32[iy][ix] - flux->data.F32[iy][ix + 1])); - } - } - - // measure max dfdy - grady = 0; - for (int iy = 0; iy < flux->numRows - 1; iy++) { - for (int ix = 0; ix < flux->numCols; ix++) { - grady = PS_MAX (grady, fabs(flux->data.F32[iy][ix] - flux->data.F32[iy + 1][ix])); - } - } - psLogMsg ("psphot.defects", 4, "gradx: %f, grady: %f\n", gradx, grady); - psFree (flux); - psFree (model); - - // search for pixels with (dfdy > grady) or (dfdx > gradx) allowing for Nsigma outliers - Nzap = 0; - for (int i = 0; i < sources->n; i++) { - psSource *source = sources->data[i]; - psImage *image = source->pixels; - psImage *mask = source->mask; - psImage *zap = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8); - psImageInit (zap, 0); - - // get sky-subtracted peak (skip sources with insufficient gradient) - sky = source->moments->Sky; - peak = source->moments->Peak - sky; - if (peak < 0) continue; - if (peak < peak*PS_MIN(gradx, grady) + NSIGMA*sqrt(2*sky + peak)) continue; - psTrace (".psphot.defects", 4, "peak: %f, SN: %f\n", peak, source->moments->SN); - - // search for (dfdx) - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols - 1; ix++) { - N0 = image->data.F32[iy][ix]; - N1 = image->data.F32[iy][ix + 1]; - gradlim = fabs(N1 - N0) - NSIGMA*(sqrt(N1 + N0)); - // gradient - N*sigma - gradlim = gradlim / peak; - - // if gradient is too large, zap the high pixel - if (gradlim < gradx) continue; - if (N0 > N1) { - zap->data.U8[iy][ix] = 1; - } else { - zap->data.U8[iy][ix + 1] = 1; - } - } - } - - // search for (dfdy) - for (int iy = 0; iy < image->numRows - 1; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - N0 = image->data.F32[iy][ix]; - N1 = image->data.F32[iy + 1][ix]; - gradlim = fabs(N1 - N0) - NSIGMA*(sqrt(N1 + N0)); - gradlim = gradlim / peak; - - // if gradient is too large, zap the high pixel - // note the inefficiency: we are re-testing known zapped entries - if (gradlim < grady) continue; - if (N0 > N1) { - zap->data.U8[iy][ix] = 1; - } else { - zap->data.U8[iy + 1][ix] = 1; - } - } - } - - // mask & zero the zapped pixels - for (int iy = 0; iy < image->numRows - 1; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (!zap->data.U8[iy][ix]) continue; - x = ix + image->col0; - y = iy + image->row0; - zapmask->data.F32[y][x] = 1; - Nzap ++; - } - } - psFree (zap); - - } - psLogMsg ("psphot", 3, "zap: %f sec, %d pixels\n", psTimerMark ("psphot"), Nzap); - - return (true); -} Index: unk/psphot/src/find_defects_new.c =================================================================== --- /trunk/psphot/src/find_defects_new.c (revision 4640) +++ (revision ) @@ -1,26 +1,0 @@ -# include "psphot.h" - -bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf) -{ - - psTimerStart ("psphot"); - - // for each source, search for all row peaks - // measure - - // 3pt second moment (assume sky?) - - // S0 = x[-1] + x[0] + x[1] - // S0 = sum (f - sky) - // S1 = sum (f - sky)*(x - x[0]) - // S2 = sum (f - sky)*(x - x[0])^2 - - // S1 = df[1] - df[-1] - // S2 = df[1]^2 + df[-1]^2 - // S0 = df[0] + df[1] + df[-1] - - // x1 = S1/S0 - // sx = sqrt (S2/S0 - x1^2) - - return (true); -} Index: /trunk/psphot/src/load_args.c =================================================================== --- /trunk/psphot/src/load_args.c (revision 4640) +++ /trunk/psphot/src/load_args.c (revision 4641) @@ -13,5 +13,5 @@ char *mask = NULL; - if ((N = get_argument (*argc, argv, "--mask"))) { + if ((N = get_argument (*argc, argv, "-mask"))) { remove_argument (N, argc, argv); mask = psStringCopy (argv[N]); @@ -20,5 +20,5 @@ char *noise = NULL; - if ((N = get_argument (*argc, argv, "--noise"))) { + if ((N = get_argument (*argc, argv, "-noise"))) { remove_argument (N, argc, argv); noise = psStringCopy (argv[N]); @@ -26,8 +26,17 @@ } + char *diff = NULL; + if ((N = get_argument (*argc, argv, "-diff"))) { + remove_argument (N, argc, argv); + diff = psStringCopy (argv[N]); + remove_argument (N, argc, argv); + } + if (*argc != 4) usage (); + // load config information psMetadata *config = psMetadataParseConfig (NULL, &Nfail, argv[3], FALSE); + // identify input image & optional noise & mask images psMetadataAddStr (config, PS_LIST_HEAD, "INPUT", "", argv[1]); psMetadataAddStr (config, PS_LIST_HEAD, "OUTPUT", "", argv[2]); @@ -38,4 +47,7 @@ psMetadataAddStr (config, PS_LIST_HEAD, "NOISE", "", noise); } + if (diff != NULL) { + psMetadataAddStr (config, PS_LIST_HEAD, "DIFF_IMAGE", "", diff); + } return (config); } @@ -43,8 +55,9 @@ int usage () { - fprintf (stderr, "USAGE: psphot (image.fits) (output.fits) (config)\n"); + fprintf (stderr, "USAGE: psphot (image) (output) (config)\n"); fprintf (stderr, "options: \n"); - fprintf (stderr, " --mask (filename)\n"); - fprintf (stderr, " --noise (filename)\n"); + fprintf (stderr, " -mask (filename)\n"); + fprintf (stderr, " -noise (filename)\n"); + fprintf (stderr, " -diff (filename)\n"); exit (2); } Index: unk/psphot/src/mark_psf_sources.c =================================================================== --- /trunk/psphot/src/mark_psf_sources.c (revision 4640) +++ (revision ) @@ -1,57 +1,0 @@ -# include "psphot.h" - -// identify objects consistent with PSF shape/magnitude distribution -// we expect dparams[4],dparams[5] to have a scatter of: -// sigma_x / (S/N) * sqrt(2) -// 1 / (params[4],params[5])*(S/N) - -// any objects which is consistent with the PSF should have -// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag) -// this should include a minimum dsx buffer for the brighter objects -// saturated stars should fall outside (but are already IDed) -// galaxies should be larger, cosmic rays smaller, but need to test? - -// need to demote PSFSTARs which don't meet criterion - -bool mark_psf_sources (psArray *sources, psMetadata *config) -{ - bool status = false; - float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); - float snFaint = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM"); - float dSX, dSY, SX, SY, SN; - float nSx, nSy; - - for (int i = 0; i < sources->n; i++) { - psSource *source = sources->data[i]; - if (source->type == PS_SOURCE_PSFSTAR) continue; - if (source->type == PS_SOURCE_SATSTAR) continue; - if (source->modelPSF == NULL) continue; - - SN = source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1]; - SX = source->modelPSF->params->data.F32[4]; - SY = source->modelPSF->params->data.F32[5]; - dSX = source->modelPSF->dparams->data.F32[4]; - dSY = source->modelPSF->dparams->data.F32[5]; - - nSx = dSX * SX * SN; - nSy = dSY * SY * SN; - - // assign PS_SOURCE_GOODSTAR to bright objects within PSF region of dparams[] - if ((fabs(nSx) < shapeNsigma) && (fabs(nSy) < shapeNsigma)) { - if (SN > snFaint) { - source->type = PS_SOURCE_GOODSTAR; - } else { - source->type = PS_SOURCE_FAINTSTAR; - } - } - - if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) { - source->type = PS_SOURCE_GALAXY; - } - if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) { - source->type = PS_SOURCE_DEFECT; - } - - } - return (true); -} Index: /trunk/psphot/src/onesource.c =================================================================== --- /trunk/psphot/src/onesource.c (revision 4640) +++ /trunk/psphot/src/onesource.c (revision 4641) @@ -14,5 +14,4 @@ float MRAD = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS"); float RADIUS = psMetadataLookupF32 (&status, config, "FIT_RADIUS"); - float SOFT = psMetadataLookupF32 (&status, config, "SOFT_RADIUS"); psSource *source = pmSourceAlloc(); Index: /trunk/psphot/src/output.c =================================================================== --- /trunk/psphot/src/output.c (revision 4640) +++ /trunk/psphot/src/output.c (revision 4641) @@ -1,27 +1,61 @@ # include "psphot.h" -// output functions -// we have several fixed modes (sx, obj, cmp) - -// we require a fixed number of chars already allocated -init_line (char *line, int *Nline) { - *Nline = 0; -} - -add_to_line (char *line, int *Nline, char *format, ...) { - - va_list ap; - - va_start (ap, format); - Nchar = vsprintf (&line[*Nline], format, ap); - *Nline += Nchar; +// output functions: we have several fixed modes (sx, obj, cmp) + +void output (psImageData *imdata, psMetadata *config, psArray *sources) { + + char *outputMode = psMetadataLookupPtr (&status, config, "OUTPUT_MODE"); + + if (outputMode == NULL) { + WriteSourcesText (imdata, config, sources); + exit (0); + } + + if (!strcasecmp (outputMode, "OBJ")) { + WriteSourcesOBJ (imdata, config, sources); + exit (0); + } + + if (!strcasecmp (outputMode, "SX")) { + WriteSourcesSX (imdata, config, sources); + exit (0); + } + + if (!strcasecmp (outputMode, "CMP")) { + WriteSourcesCMP (imdata, config, sources); + exit (0); + } + + if (!strcasecmp (outputMode, "CMF")) { + WriteSourcesCMF (imdata, config, sources); + exit (0); + } + + psAbort ("psphot", "unknown output mode %s", outputMode); +} + +bool WriteSourceText (psImageData *imdata, psMetadata *config, psArray *sources) { + + char *filename = psMetadataLookupPtr (&status, config, "OUTPUT"); + + DumpImage (imdata->image, filename); + + // get these names from config? + DumpModelPSF (sources, "psfsources.dat"); + DumpModelFLT (sources, "fltsources.dat"); + DumpModelNULL (sources, "nullsources.dat"); + DumpMoments (sources, "moments.dat"); + + exit (0); } // dophot-style output list with fixed line width -bool WriteSourcesOBJ (psArray *sources, char *filename) { +bool WriteSourcesOBJ (psImageData *imdata, psMetadata *config, psArray *sources) { # define NCHAR 104 char line[NCHAR]; + + char *filename = psMetadataLookupPtr (&status, config, "OUTPUT"); FILE *f = fopen (filename, "w"); @@ -63,9 +97,11 @@ // elixir/sextractor-style output list with fixed line width -bool WriteSourcesSX (psArray *sources, char *filename) { +bool WriteSourcesSX (psImageData *imdata, psMetadata *config, psArray *sources) { # define NCHAR 104 char line[NCHAR]; + + char *filename = psMetadataLookupPtr (&status, config, "OUTPUT"); FILE *f = fopen (filename, "w"); @@ -108,5 +144,5 @@ // elixir-style pseudo FITS table (header + ascii list) -bool WriteSourcesCMP (psMetadata *header, psArray *sources, char *filename) { +bool WriteSourcesCMP (psImageData *imdata, psMetadata *config, psArray *sources) { @@ -114,4 +150,6 @@ char line[NCHAR]; + + char *filename = psMetadataLookupPtr (&status, config, "OUTPUT"); FILE *f = fopen (filename, "w"); @@ -121,5 +159,5 @@ } - // write header to file + // write imdata->header to file // write sources with models first @@ -155,5 +193,25 @@ // elixir-style FITS table output (header + table in 1st extension) -bool WriteSourcesCMF (psArray *sources, char *filename) { - // ??? -} +bool WriteSourcesCMF (psImageData *imdata, psMetadata *config, psArray *sources) { + + char *filename = psMetadataLookupPtr (&status, config, "OUTPUT"); + + // write imdata->header to file + + // write FITS table to file (use autocode tools) +} + +// we require a fixed number of chars already allocated +init_line (char *line, int *Nline) { + *Nline = 0; +} + +add_to_line (char *line, int *Nline, char *format, ...) { + + va_list ap; + + va_start (ap, format); + Nchar = vsprintf (&line[*Nline], format, ap); + *Nline += Nchar; +} + Index: unk/psphot/src/pixregion.c =================================================================== --- /trunk/psphot/src/pixregion.c (revision 4640) +++ (revision ) @@ -1,14 +1,0 @@ -# include "psphot.h" - -psRegion *define_valid_pixels (psImage *image, psMetadata *config) -{ - - bool status; - float XBORDER = psMetadataLookupF32 (&status, config, "XBORDER"); - float YBORDER = psMetadataLookupF32 (&status, config, "YBORDER"); - - psRegion *pixregion = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER); - pixregion = psRegionForImage (pixregion, image, pixregion); - - return (pixregion); -} Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 4640) +++ /trunk/psphot/src/psphot.c (revision 4641) @@ -34,7 +34,4 @@ psf = choose_psf_model (config, sources, sky); - // identify defects based on 1-pixel wide x and y direction 2nd moments - // find_defects (zap, sources, config, psf); - // test PSF on all except SATURATE and DEFECT (artifacts) apply_psf_model (imdata, config, sources, psf, sky); @@ -43,11 +40,6 @@ fit_galaxies (imdata, config, sources, sky); - // subtract_galaxies (sources, config); - - DumpImage (imdata->image, argv[2]); - DumpModelPSF (sources, "psfsources.dat"); - DumpModelFLT (sources, "fltsources.dat"); - DumpModelNULL (sources, "nullsources.dat"); - DumpMoments (sources, "moments.dat"); + // write out data in appropriate format + output (imdata, config, sources); exit (0); Index: /trunk/psphot/src/setup.c =================================================================== --- /trunk/psphot/src/setup.c (revision 4640) +++ /trunk/psphot/src/setup.c (revision 4641) @@ -44,4 +44,6 @@ } else { + // XXX EAM it looks like psBinaryOp is broken for: image op scalar + // this loop performs the operation by hand noise = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32); for (int iy = 0; iy < image->numRows; iy++) { @@ -51,5 +53,5 @@ } - // XXX EAM it looks like psBinaryOp is broken for (image op scalar) + // this uses psBinaryOp # if (0) psScalar *value; @@ -65,10 +67,11 @@ } + // load the mask if specified, otherwise construct char *maskName = psMetadataLookupPtr (&status, config, "MASK"); if (status == true) { file = psFitsAlloc (maskName); mask = psFitsReadImage (NULL, file, region, 0); - // require U8 psFree (file); + // XXX require / convert mask to psU8? } else { mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8); @@ -92,8 +95,7 @@ } } - DumpImage (mask, "mask.fits"); - psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot")); + // save the image data & return it psImageData *imdata = psAlloc(sizeof(psImageData)); imdata->image = image; Index: unk/psphot/src/subtract_galaxies.c =================================================================== --- /trunk/psphot/src/subtract_galaxies.c (revision 4640) +++ (revision ) @@ -1,23 +1,0 @@ -# include "psphot.h" - -// fit selected galaxy model (GAUSS) to all bright objects of type GALAXY - -bool subtract_galaxies (psArray *sources, psMetadata *config) -{ - for (int i = 0; i < sources->n; i++) { - psSource *source = sources->data[i]; - if (source->type != PS_SOURCE_GALAXY) continue; - - psModel *model = source->modelFLT; - if (model == NULL) continue; - - psImage *pixels = source->pixels; - - // subtract object, leave local sky - psF64 sky = model->params->data.F32[0]; - model->params->data.F32[0] = 0; - pmSourceSubModel (pixels, source->mask, model, false); - model->params->data.F32[0] = sky; - } - return (true); -} Index: unk/psphot/src/subtract_psf_sources.c =================================================================== --- /trunk/psphot/src/subtract_psf_sources.c (revision 4640) +++ (revision ) @@ -1,42 +1,0 @@ -# include "psphot.h" - -// subtract model fits from image (only PSFSTAR && BRIGHTSTAR) -// need to control application of this with some thresholding - -bool subtract_psf_sources (psArray *sources) -{ - - sources = psArraySort (sources, by_SN); - - for (int i = 0; i < sources->n; i++) { - float sky; - psSource *source = sources->data[i]; - - if (source->type == PS_SOURCE_SATSTAR) { - if (source->modelPSF == NULL) { - fprintf (stderr, "missing model for %f, %f\n", source->moments->x, source->moments->y); - continue; - } - fprintf (stderr, "subtracting model for %f, %f\n", source->moments->x, source->moments->y); - } - - - bool valid = false; - valid |= (source->type == PS_SOURCE_SATSTAR); - valid |= (source->type == PS_SOURCE_PSFSTAR); - valid |= (source->type == PS_SOURCE_BRIGHTSTAR); - if (!valid) continue; - - psModel *model = source->modelPSF; - if (model == NULL) continue; - - psImage *pixels = source->pixels; - - // subtract object, leave local sky - sky = model->params->data.F32[0]; - model->params->data.F32[0] = 0; - pmSourceSubModel (pixels, source->mask, model, false); - model->params->data.F32[0] = sky; - } - return(true); -}