Index: /trunk/archive/scripts/phase2.c =================================================================== --- /trunk/archive/scripts/phase2.c (revision 2519) +++ /trunk/archive/scripts/phase2.c (revision 2519) @@ -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/scripts/phase2.pl =================================================================== --- /trunk/archive/scripts/phase2.pl (revision 2519) +++ /trunk/archive/scripts/phase2.pl (revision 2519) @@ -0,0 +1,93 @@ +#!/usr/bin/env perl + +sub usage { + print STDERR "USAGE: phase2.pl\n"; + exit 2; +} + +if (@ARGV != 4) { &usage; } + +# phase2.pl has extensive intelligence about data sources + +# input: filename.fits +# output: flattened image: filename.f.fits? filename.flt? (chosen externally?) +# output: extracted stars: filename.p.tbl +# output: constructed kernel: filename.k.fits + +# determine camera from image or externally? + +# load configuration information: + +$config = psMetadataReadFile (NULL, $configfile); + +# check image: +# is it a FITS image? +# find camera +# find mode: +# one cell per extension +# one chip per extension +# one FPA per extension + +@chips = getChips ($input); +foreach $chip () { + @cells = getCells ($input, $chip); +} + + + +# construct kernel + +# load bias +# load dark +# apply bias and dark + +# all possible detrend operations +# - bias +# - dark +# - trim +# - mask bad pixels (astrometry independent) +# - mask artifacts (astrometry dependent) +# - flatten +# - fringe correct + +# convolution kernel is applied to: +# - flat +# - fringe +# - bias? + +# + + +$out = pmSubtractBias ($in, $fit); + + +##### example: single chip, single cell, single readout: + +# input filename: $input + +# bias recipe: +# bias image: filename +# overscan: header BIASSE + +# bias subtraction options: +# - use a bias image? +# - + + + +$image = psImageReadSection (NULL, 0, 0, 0, 0, 0, $extname, -1, $input); +$metadata = psMetadataReadHeader (NULL, $extname, -1, $input); + +$metaitem = psMetadataLookup ($metadata, "BIASSEC"); +$overscan = psImageSubsection ($image, $metaitem->data.void); + +psListAdd ($overlist, $overscan, PS_LIST_HEAD); + +$bias = psImageReadSection (NULL, 0, 0, 0, 0, 0, $extname, -1, $biasname); + +$image_b = pmSubtractBias ($image, NULL, $overlist, PM_OVERSCAN_COLUMNS, $stats, 0, PM_FIT_NONE, $bias); + + + +####################### + Index: /trunk/archive/scripts/phase4.c =================================================================== --- /trunk/archive/scripts/phase4.c (revision 2519) +++ /trunk/archive/scripts/phase4.c (revision 2519) @@ -0,0 +1,212 @@ +#include +#include + +// Configuration for Phase 4 +typedef struct { + // Input parameters + int nInputs; // Number of input files + const char *camera; // Name of the camera configuration file + const psList *inNames; // Input file names + unsigned long skyId; // Sky tangent plane identification number + // Output parameters + const char *combined; // Combined image filename + const char *subtracted; // Subtracted image filename + const char *variables; // Output variable candidates filename + const char *recipeFile; // Recipe filename +} phase4Config; + + +// A tangent plane (for combined image and static sky) +// Assume we're dealing with a rectangle only for now. +typedef struct { + long int id; // Identification number --- key in MDDB + double ra0, dec0; // Tangent point --- from MDDB + double sizeRA, sizeDec; // Size in RA and Dec --- from MDDB + double scale; // Pixel scale --- from MDDB + char *name; // File name --- from MDDB or image server + psImage *pixels; // The pixels that comprise the tangent plane + psImage *mask; // Mask image + psSomeCollectionOfObjects *objects; // Objects in the image +} psTangentPlane; + +void main (int argc, char **argv) +{ + // Bury the command line parsing into some function that just works + phase4Config *config = parseCommandLine(argc, argv); // Parse the command line + + // Read the recipe metadata + psMetadata *recipe = psMetadataParseConfig(NULL, config->recipeFile, false); + + // Don't particularly care about the camera and site details --- Phase 4 doesn't care where + // the data comes from, but it does need the following for each image: + // * Gain and read noise for noise model --- from the camera configuration + // * Background value for noise model (assume constant over the region of interest) --- from image MD + // * Photometric calibration for relative and absolute scalings --- from image MD + // * Astrometric calibration to map to the sky --- from file specified in image MD + // We need to be told which camera; assume all data is from the same camera (easy to generalise + // if this is not true). + + // Read the camera configuration + psMetadata *camera = psMetadataParseConfig(NULL, config->camera, false); + // What are the keywords for gain, readnoise, photometric calibration, background? + // This basically treats the "camera" configuration as a translation table. + const char *gainHeader = psMetadataLookup(camera, "GAIN"); + const char *readnoiseHeader = psMetadataLookup(camera, "READNOISE"); + const char *photcalHeader = psMetadataLookup(camera, "PHOTCAL"); + const char *backgroundHeader = psMetadataLookup(camera, "BACKGRND"); + const char *astromHeader = psMetadataLookup(camera, "ASTROM"); + + // Initialise the various arrays/vectors we need + psArray *headers = psArrayAlloc(config->nInputs); // The image headers + psArray *inputs = psArrayAlloc(config->nInputs); // The images + psArray *astrometry = psArrayAlloc(config->nInputs); // The astrometry data for each image + psVector *gain = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Gains + psVector *readnoise = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Read noises + psVector *photcal = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Photometric calibrations + psVector *background = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Background flux values + + // This function queries the MDDB to get sizeRA,sizeDec,ra0,dec0,scale,name for the tangent plane of + // interest + psTangentPlane *staticSky = psTangentPlaneFromDB(config->sky); + + // Read in each of the inputs + for (psListSetIterator(config->inNames, PS_LIST_HEAD), // Set the iterator and get the first one + char *inName = psListGetNext(input), + int i = 0; + inName != NULL; // Keep going until nothing left in the list + inName = psListGetNext(input), // Get the next one on the list + i++) { + + // Read the header + psMetadata *header = psFITSReadHeader(inName); + // Put it in the array of headers for future reference + headers->data[i] = header; + + // Check the validity of the given header based on expectactions for this camera + status = pmCameraValidateHeaders(headers, camera); + // Need to act if status is bad --- generate a warning and ignore this image + + // Get the gain, readnoise, photometric calibration, background + // I think this is not strictly correct --- doesn't psMetadataLookup return a void* ? + gain->data.F32[i] = (psF32) psMetadataLookup(header, gainHeader); + readnoise->data.F32[i] = (psF32) psMetadataLookup(header, readnoiseHeader); + photcal->data.F32[i] = (psF32) psMetadataLookup(header, photcalHeader); + background->data.F32[i] = (psF32) psMetadataLookup(header, backgroundHeader); + + // Get the name of the astrometry file (XML) and generate the FPA structure from it (no pixels) + const char *astromFile = psMetadataLookup(header, astromHeader); + psFPA *fpa = psFPAFromAstrometry(astromFile); + astrometry->data[i] = fpa; + + // Iterate over the FPA structure + for (int chipnum = 0, int nChips = fpa->chips->n; chipnum < nChips; chipnum++) { + psChip *chip = fpa->chips->data[chipnum]; // The chip of interest + for (int cellnum = 0, int nCells = chip->cells->n; cellnum < nCells; cellnum++) { + psCell *cell = chip->cells->data[cellnum]; // The cell of interest + + // QUESTION: Do we want to bother about the cells that were used for fast guiding? + // Assume that we care --- may as well. + + // QUESTION: Are multiple readouts stacked in phase 4 or earlier? + // Assume that they should be stacked here. + + // QUESTION: How do we know the names of the chip, cell, readout in the FITS file? + // Assume that each of psFPA, psChip, psCell, psReadout contains a "const char *source" + // which contains some sort of URI, e.g.: + // FPA: fpa12345_12.fits --- This particular FPA contains only a single OTA + // Chip: fpa12345_12.fits --- The OTA is entirely contained in the FITS file + // Cell: fpa12345_12.fits|cell34 --- The cell has extension name "cell34" + // Readout: fpa12345_12.fits|cell34|56 --- The readout is the 56th slice in the 3rd axis + // Assume further that, given the "source", there exists a module that loads the pixels. + + + // If there are multiple readouts, then stack them + if (cell->readouts->n > 1) { + psList *readouts = psListAlloc(NULL); // List of readouts for stacking + + // Read the images in + for (int readoutNum = 0; readoutNum < cell->readouts->n; readoutNum++) { + (void)psReadoutRead(cell->readouts->data[readoutNum]); // Read the readout into memory + psListAdd(readouts, PS_LIST_TAIL, cell->readouts->data[readoutNum]); // Put on list + } + + // QUESTION: Do we want to bother about combining these "properly", or just sum them? + // Assume that we just sum them, which saves measuring the zeros and scales. CR rejection + // can be done when we combine multiple images (consecutive or from multiple telescopes). + + // Perhaps we want a pmReadoutStack module, which would convert an array of readouts from + // fast guiding to a single pseudo-readout. + + // Just do a simple mean + psStats *combineStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics for combination + pmCombineParams *combine = pmCombineParamsAlloc(combineStats, 0, 0.0, 0.0, 0); // How to + // combine + psImage *stack = pmReadoutCombine(NULL, readouts, combine, 0.0, 1.0, false, 1.0, 0.0); + // And now multiply by the number to get the sum + (void)psImageBinaryOp(stack, stack, "*", psScalarAlloc(cell->readouts->n, PS_TYPE_U16)); + + // Probably also have to multiply by some exposure time ratio to scale with the rest of + // the focal plane --- neglected here. + + // Clean up the temporary steps + psFree(readouts); + psFree(combine); + psFree(combineStats); + + // Stuff the stack back into the cell --- probably should play with the metadata, but I'm + // going to ignore that for now. + + // Clean out the old structure + for (int readoutNum = 0; readoutNum < cell->readouts->n; readoutNum++) { + psFree(cell->readouts->data[readoutNum]); + } + // Resize and stuff the stack in. + psArrayRealloc(cell->readouts, 1); + cell->readouts->data[0] = stack; + } else { + // Otherwise, just read the cell in + cell->readouts->data[0] = psReadoutRead(cell->readouts->data[0]); + } // Reading in readouts + } // Cells + } // Chips + } // FPAs + + // Now all the data is in memory. Combine the multiple images + psTangentPlane *combined = pmFPACombine(NULL, staticSky, inputs, gain, readnoise, photcal, background); + + // Write the combined image out + psTangentPlaneWriteFITS(combined, config->combined); + + // Run the object detection/measurement --- this is not well thought out yet + (void)psDetectAndMeasureObjects(combined, otherParameters); + + // Load the static sky into the structure already defined + (void)psTangentPlaneRead(staticSky); + + // Get the seeing ratio --- involves matching stars and measuring the mean ratio of FWHM + float seeingRatio = pmSeeingRatio(combined, staticSky, otherParameters); + + // Do the subtraction --- involves measuring the convolution kernel, applying the kernel and subtracting + psTangentPlane *subtracted = NULL; + if (seeingRatio < 1.0) { + subtracted = pmImageSubtraction(NULL, combined, staticSky, otherParameters); + } else { + subtracted = pmImageSubtraction(NULL, staticSky, combined, otherParameters); + } + + // Write out the subtracted image + psTangentPlaneWriteFITS(subtraction, config->subtracted); + + // Find and measure objects on the subtracted image + (void)psDetectAndMeasureObjects(subtracted, otherParameters); + + // Filter the list of candidate variable sources + psList *variables = pmFilterVariables(subtracted, combined, staticSky, otherParameters); + + // Output the list of variable sources + (void)pmOutputVariables(variables, config->variables); + + // All done. +} + + Index: /trunk/archive/scripts/pmFPAio.c =================================================================== --- /trunk/archive/scripts/pmFPAio.c (revision 2519) +++ /trunk/archive/scripts/pmFPAio.c (revision 2519) @@ -0,0 +1,20 @@ +# include + +main (int argc, char **argv) { + + psArray headers; + psMetadata *base, *site, *camera, *recipe; + + pmConfigLoadBasic (&site, &base); + + fits = psFitsAlloc (input); + header = psFitsReadHeader (fits); + pmConfigLoadCamera (&camera, site, header); + + pmConfigLoadRecipe (&recipe, camera, "PHASE2"); + + // we've loaded the PHU, determined the camera (and recipe). + + fpa = psFPAAlloc (); + + // is this an image?