Index: /trunk/archive/modules/src/phase2.c =================================================================== --- /trunk/archive/modules/src/phase2.c (revision 2069) +++ /trunk/archive/modules/src/phase2.c (revision 2069) @@ -0,0 +1,300 @@ +# include + +// ignoring things like allocation, error tests, etc + +main (int argc, char **argv) { + + psArray headers; + psMetadata *base, *site, *camera, *recipe; + + // the set of code below (up to END CONFIG) is generic boilerplate which can be used for most analysis stages + // the specific issues are defined by the recipefile loaded for this analysis + // I've written this keeping the external metadata separated (base, site, camera, recipe), + // but we could put all of this in a single metdata collection, called eg config + + // USAGE: phase2 (input) (outroot) [-camera (camera)] [-recipe (recipe)] + input = argv[1]; + outroot = argv[2]; + + // define command-line options + cameraname = NULL; + cameraname = getargs ("-camera"); // pseudo-code, not a valid API + recipefile = NULL; + recipefile = getargs ("-recipe"); // pseudo-code, not a valid API + + // base file needs to be hardwired, from environment, or based on program name + base = psMetadataParseConfig (NULL, basefile, FALSE); + + // load site-specific information (such as database servers, etc) + sitefile = psMetadataLookup (base, "SITE"); + site = psMetadataParseConfig (NULL, sitefile, FALSE); + + // load only the headers from the given file input a list of headers + headers = psFITSReadHeaderSet (input); + + // identify camera (may be specified on command line) + if (cameraname == NULL) { + cameraname = pmCameraFromHeader (headers.data[0]); + } + camerafile = psMetadataLookup (site, "CAMERA.%s", cameraname); + + // load the camera details for the specific camera + camera = psMetadataParseConfig (NULL, camerafile, FALSE); + + // recipefile is defined in the camera metadata, if not specified on command line + if (recipefile == NULL) { + recipefile = psMetadataLookup (camera, "RECIPE.PHASE2"); + recipefile = DetrendDatabaseLookup ("recipe", camera, header); // an alternative + // make this a database lookup & a function of date + } + recipe = psMetadataParseConfig (NULL, recipefile, FALSE); + + // END CONFIG loading code + + // check the validity of the given set of headers based on expectactions for this camera + status = pmCameraValidateHeaders (headers, camera); + + // construct a complete, empty (ie, no pixel data) FPA structure for this camera / header set + // we need to construct the output container. is it big enough to fit in memory? assume so? + FPA = pmFPAfromHeader (headers, camera); + + // I'm pretending the extname value for each cell in part of psCell + // Do we need to construct a list of the extname values associated with each readout? + // Is EXTNAME invarient? (probably not; does not always exist) + + for (i = 0; i < FPA[0].chips.n; i++) { + chip = FPA[0].chips.data[i]; + + for (j = 0; j < chip[0].cells.n; j++) { + cell = chip[0].cells.data[j]; + + bias = load_bias_image (recipe, cell[0].metadata, cell[0].extname); + dark = load_dark_image (recipe, cell[0].metadata, cell[0].extname); + flat = load_flat_image (recipe, cell[0].metadata, cell[0].extname); + badpix = load_badpix_image (recipe, cell[0].metadata, cell[0].extname); + fringe = load_fringe_image (recipe, cell[0].metadata, cell[0].extname); + + for (k = 0; k < cell[0].readouts.n; k++) { + readout = cell[0].readouts[k]; + + // load the specific readout pixels. put this in a function? + readout[0].image = psImageReadSection (NULL, 0, 0, 0, 0, 0, cell[0].extname, k, input); + + kernel = make_kernel (); + mask = make_mask (); + + bias_subtract_readout (readout, bias, recipe); + dark_subtract_readout (readout, dark, recipe); + nonlinear_correct_readout (readout, mask, recipe); + + trim_readout (readout, mask, recipe); + mask_readout (readout, mask, badpix, kernel, recipe); + + // where do we convolve the flat with the kernel? + pmFlatField (readout, mask, flat); + + defringe_readout (readout, mask, fringe, recipe); + } + } + } +} + +// bias subtract the given readout according to the recipe (bias image may be NULL) +// no mask needed yet +bias_subtract_readout (psReadout *readout, psImage *bias, psMetadata *recipe) { + + char *biassec; + psImage *overscan; + psList *overscanlist; + pmOverscanAxis direction; + psRectangle biassec; + + // determine overscan option + biassec = pmRecipeDefineSection (recipe, readout[0].metadata, "BIAS.OVERSCAN"); + overscan = psImageSubsection (readout[0].image, biassec); + + psListAdd (overscanlist, overscan, PS_LIST_HEAD); // list of overscans (allocate list first) + + // the direction is implied by noting the long dimension of the biassec + direction = DirectionFromOverscan (biassec); // weak, but do we want a function to deal with this or add it inline? + + // the type of statistic is defined by the recipe: BIAS.OVERSCAN.STATS + statname = psMetadataLookup (recipe, "BIAS.OVERSCAN.STATS"); // string defining overscan stats + stats = DefineStats (statname); + + // identify the type of fit requested and relevant parameters + // the type of fit is defined by the recipe: BIAS.OVERSCAN.STATS + fitname = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT"); // string defining overscan stats + + if (!strcmp (fitname, "NONE")) { + fitType = PM_FIT_NONE; + fitSpec = NULL; + nBin = 0; + } + + if (!strcmp (fitname, "SPLINE")) { + fitType = PM_FIT_SPLINE; + nBin = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.NBIN"); + nPts = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.NPTS"); + order = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.ORDER"); // only 1 or 3 is allowed for 'order' + fitSpec = psSpline1DAlloc (nPts, 3, 0, Npix); // Npix is either NAXIS1 or NAXIS2 (Nx or Ny); + } + + if (!strcmp (fitname, "POLYNOMIAL")) { + fitType = PM_FIT_POLYNOMIAL; + nBin = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.NBIN"); + nPts = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.NPTS"); + order = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.ORDER"); + fitSpec = psPolynomial1DAlloc (nPts, PS_POLYNOMIAL_CHEB); + } + + readout[0].image = pmSubtractBias (readout[0].image, fitSpec, overscanlist, direction, stats, nBin, PM_FIT_NONE, bias); + + // add these to the readout[0].metadata? + // or to the chip / cell metadata? + +} + +// subtract dark image using the appropriate scaling factor for this exposure time +// no mask needed yet +dark_subtract_readout (psReadout *readout, psImage *dark, psMetadata *recipe) { + + char *biassec; + psImage *overscan; + psList *overscanlist; + pmOverscanAxis direction; + + if (dark == NULL) { + return; + } + + // need to get a lookup table from the database, apply exptime for readout and dark + factor_source = psMetadataLookup (recipe, "DARK.FACTOR"); // use the recipe to determine the dark scaling + factor = getScalingFactor (readout, dark, factor_source); + + // out = in - dark * factor + tmpdark = psBinaryOp (NULL, dark, "*", factor); + readout[0].image = psBinaryOp (readout[0].image, readout[0].image, "-", tmpdark); + + psFree (tmpdark); +} + +// this corrects the non-linear response, but probably should use data determined at the cell level +// we need to set a mask value for pixels which have were out of range (TBD) +nonlinearity_correct_readout (psReadout *readout, psReadout *mask, psMetadata *recipe) { + + char *method; + psPolynomial1D correction; + psVector inFlux, outFlux; + + method = psMetadataLookup (recipe, "NONLINEAR.METHOD"); + + // use a polynomial; get coeffs from recipe + if (!strcmp (method, "POLYNOMIAL")) { + correction.n = psMetadataLookup (recipe, "NONLINEAR.ORDER"); + correction.coeffs = psMetadataLookup (recipe, "NONLINEAR.COEFFS"); + pmNonLinearityPolynomial (readout, correction); + } + + if (!strcmp (method, "LOOKUP")) { + lookup_table = DetrendDatabaseLookup (time, cell, camera, etc?); + inFlux = psFITSTableRead (lookup_table, "IN_FLUX"); + outFlux = psFITSTableRead (lookup_table, "OUT_FLUX"); + pmNonLinearityLookup (readout, inFlux, outFlux); + } +} + +// trim the image and mask based on the requested region +trim_readout (psReadout *readout, psReadout *mask, psMetadata *recipe) { + + datasec = pmRecipeDefineSection (recipe, readout[0].metadata, "TRIM.DATASEC"); + + psImageTrim (readout[0].image, datasec); + psImageTrim (mask[0].image, datasec); + +} + +// apply the bad pixel mask to the image mask +mask_readout (psReadout *readout, psImage *mask, psImage *badpix, psKernel *kernel, psMetadata *recipe) { + +} + +psImage *load_bias_image (psMetadata *recipe, psMetadata *header, char *extname) { + + // load in the complete corresponding bias image (this should be done once per cell) + // we are implying that the bias image is guaranteed to have the corresponding EXTNAME + bias_source = psMetadataLookup (recipe, "BIAS.IMAGE"); + biasname = NULL; + if (!strncmp (bias_source, "FILE:", 5)) { + biasname = &bias_source[5]; + } + if (!strcmp (bias_source, "DB:BEST")) { + biasname = DetrendDatabaseLookup ("best", header); // these APIs need to be seriously fleshed out!! + } + if (!strcmp (bias_source, "DB:CLOSE")) { + biasname = DetrendDatabaseLookup ("close", header); // these APIs need to be seriously fleshed out!! + } + if (biasname != NULL) { + bias = psImageReadSection (NULL, 0, 0, 0, 0, 0, extname, 0, biasname); + } + + return (bias); +} + +psImage *load_dark_image (psMetadata *recipe, psMetadata *header, char *extname) { + + // load in the complete corresponding bias image (this should be done once per cell) + // we are implying that the bias image is guaranteed to have the corresponding EXTNAME + dark_source = psMetadataLookup (recipe, "DARK.IMAGE"); + darkname = NULL; + if (!strncmp (dark_source, "FILE:", 5)) { + darkname = &dark_source[5]; + } + if (!strcmp (dark_source, "DB:BEST")) { + darkname = DetrendDatabaseLookup ("best", header); // these APIs need to be seriously fleshed out!! + } + if (!strcmp (dark_source, "DB:CLOSE")) { + darkname = DetrendDatabaseLookup ("close", header); // these APIs need to be seriously fleshed out!! + } + if (darkname != NULL) { + dark = psImageReadSection (NULL, 0, 0, 0, 0, 0, extname, 0, darkname); + } + + return (dark); +} + +psRectangle *pmRecipeDefineSection (psMetadata *recipe, psMetadata *header, char *name) { + + char *source, *section; + psRectangle *Section; + + source = psMetadataLookup (recipe, name); + + if (!strcasecmp (source, "NONE")) { + Section = psRectangleAlloc (0, 0, 0, 0); // the NULL section + return (Section); + } + + if (!strncasecmp (source, "HEADER:", 7)) { + char *keyword; + + keyword = &source[7]; + section = psMetadataLookup (readout[0].metadata, keyword); + Section = pmSectiontoRectangle (section); + psFree (section); + psFree (source); + return (Section); + } + + if (!strncasecmp (source, "RECIPE:", 7)) { + section = &source[7]; + Section = pmSectiontoRectangle (section); + psFree (source); + return (Section); + } + + psError ("invalid section"); + return (NULL); + +} + Index: /trunk/archive/modules/src/phase2.txt =================================================================== --- /trunk/archive/modules/src/phase2.txt (revision 2069) +++ /trunk/archive/modules/src/phase2.txt (revision 2069) @@ -0,0 +1,134 @@ + +# these are examples of camera definition variables: +# skyprobe +NCELL S32 1 +NCHIP S32 1 +# FILENAME EXTNAME REGION CHIP BIASSEC +CELL.00 STR %f00.%x PHU [0,0:0,0] CHIP.00 [0,0:0,0] + +# megacam-raw +NCELL S32 72 +NCHIP S32 36 +# FILENAME EXTNAME DATASEC CHIP BIASSEC +CELL.00 STR %f.%x AMP00 [0,0:0,0] CHIP.00 BIASSEC +CELL.01 STR %f.%x AMP01 [0,0:0,0] CHIP.00 [2100,2110:0,4096] +CELL.02 STR %f.%x AMP02 [0,0:0,0] CHIP.01 [0,0:0,0] +CELL.03 STR %f.%x AMP03 [0,0:0,0] CHIP.01 [0,0:0,0] + +# megacam-splice +NCELL S32 72 +NCHIP S32 36 + +TYPE CELL FILENAME:STR EXTNAME REGION BIASSEC TRIMSEC CHIP +CELL.00 CELL %f.%x CCD00 HEADER:ASEC-00 BSEC-00 DSEC-00 CHIP.00 +CELL.01 CELL %f.%x CCD00 HEADER:ASEC-01 BSEC-01 DSEC-01 CHIP.00 +CELL.02 CELL %f.%x CCD01 HEADER:ASEC-00 BSEC-00 DSEC-00 CHIP.01 +CELL.03 CELL %f.%x CCD01 HEADER:ASEC-01 BSEC-01 DSEC-01 "chip foo" +CELL.03 CELL DEFAULT=CELL.00 FILENAME=%f.%x + +FPA: + CHIP.00 + CELL.00 CELL %f.%x CCD00 ASEC-00 BSEC-00 DSEC-00 CHIP.00 + CELL.01 CELL %f.%x CCD00 ASEC-01 BSEC-01 DSEC-01 CHIP.00 + CHIP.01 + CELL.02 CELL %f.%x CCD01 ASEC-00 BSEC-00 DSEC-00 CHIP.01 + CELL.03 CELL %f.%x CCD01 ASEC-01 BSEC-01 DSEC-01 "chip foo" + + +# gpc +NCELL S32 4 +NCHIP S32 2 +# FILENAME EXTNAME DATASEC CHIP BIASSEC +CHIP.00.CELL.00 STR %f/%f00.%x CELL00 "[0,0:0,0]" CHIP.00 RECIPE:[0,0:0,0] +CHIP.00.CELL.01 STR %f/%f00.%x CELL01 "DATASEC" CHIP.00 RECIPE:[0,0:0,0] +# + +TYPE CHIP NCELL +TYPE CELL FILENAME:STR EXTNAME REGION BIASSEC TRIMSEC CHIP + +CHIP.01 CHIP 2 +# TYPE CELL FILENAME:STR EXTNAME REGION BIASSEC TRIMSEC CHIP + CHIP.01.CELL.00 CELL %f/%f01.%x CELL02 RECIPE:[0,0:0,0] CHIP.01 [0,0:0,0] + CHIP.01.CELL.01 CELL %f/%f01.%x CELL03 RECIPE:[0,0:0,0] CHIP.01 [0,0:0,0] + +CHIP.00 CHIP + + + +f12345 +CHIP37 + + + +# cfh12k-split +NCELL S32 12 +NCHIP S32 12 +# FILENAME EXTNAME REGION CHIP BIASSEC +CELL.00 STR %f/%f00.%x PHU [0,0:0,0] CHIP.00 [0,0:0,0] +CELL.01 STR %f/%f01.%x PHU [0,0:0,0] CHIP.01 [0,0:0,0] +CELL.02 STR %f/%f02.%x PHU [0,0:0,0] CHIP.02 [0,0:0,0] + +# cfh12k-mef +NCELL S32 12 +NCHIP S32 12 +# FILENAME EXTNAME REGION CHIP BIASSEC +CELL.00 STR %f.%x CHIP00 [0,0:0,0] CHIP.00 [0,0:0,0] +CELL.01 STR %f.%x CHIP01 [0,0:0,0] CHIP.01 [0,0:0,0] +CELL.02 STR %f.%x CHIP02 [0,0:0,0] CHIP.02 [0,0:0,0] + +- REGION can be defined by a header keyword in IRAF format or by a explicit IRAF format +- what is default for NAXIS1,2 for IRAF format? +- Nreadout is always NAXIS3? + +# recipe file: +# this makes the assumption that, for a given camera, all chips & +# cells have the same recipe. this is probably a good start, but may +# not cut it in general. eg, it is already clear that for + +# recipe file must be a function of time and camera. +# + +# BIAS: +BIAS.IMAGE STR NONE +BIAS.IMAGE STR FILE:bias.fits +BIAS.IMAGE STR DB:BEST +BIAS.IMAGE STR DB:CLOSE + +BIAS.OVERSCAN STR HEADER:BIASSEC +BIAS.OVERSCAN STR RECIPE:[0,0:0,0] +BIAS.OVERSCAN STR NONE + +BIAS.OVERSCAN.STATS STR MEDIAN +BIAS.OVERSCAN.STATS STR MEAN + +BIAS.OVERSCAN.FIT STR SPLINE +BIAS.OVERSCAN.FIT.NPTS S32 5 + +BIAS.OVERSCAN.FIT STR POLYNOMIAL +BIAS.OVERSCAN.FIT.ORDER S32 3 +BIAS.OVERSCAN.FIT.NBIN S32 5 + +bias image: + none + filename + database: best + database: acceptable match + +overscan: + from header: BIASSEC + use region ['x0,x1:y0:y1'] + statistics: median + direction PM_OVERSCAN_COLUMNS + fit: POLYNOMIAL + order: 3 + +bias +dark +non-linearity +flat +mask +trim +fringe (no large-scale structure?) + +detect objects +astrometrize objects