Index: /trunk/archive/scripts/Makefile =================================================================== --- /trunk/archive/scripts/Makefile (revision 3703) +++ /trunk/archive/scripts/Makefile (revision 3703) @@ -0,0 +1,71 @@ +default: psphot +help: + @echo "USAGE: make psphot" + +CC = gcc -g -Wall +SRC = src +BIN = bin + +DESTBIN = /home/eugene/src/bin/$(ARCH) + +LPSLIB := $(shell pslib-config --libs) +IPSLIB := $(shell pslib-config --cflags) + +INCS = $(IPSLIB) +LIBS = -lpsmodule $(LPSLIB) +CFLAGS = $(INCS) +LFLAGS = $(LIBS) -static + +PSPHOT = \ +$(SRC)/psphot.$(ARCH).o + +psphot: $(BIN)/psphot.$(ARCH) + +$(BIN)/psphot.$(ARCH) : $(PSPHOT) + +# $(PSPHOT): $(INC)/psphot.h + +PSTEST = \ +$(SRC)/pstest.$(ARCH).o + +pstest: $(BIN)/pstest.$(ARCH) + +$(BIN)/pstest.$(ARCH) : $(PSTEST) + +INSTALL = psphot + +# dependancy rules for binary code ######################### +.PRECIOUS: %.$(ARCH).o +.PRECIOUS: $(BIN)/%.$(ARCH) + +%.$(ARCH).o : %.c + $(CC) $(CFLAGS) -c $< -o $@ + +$(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o + @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi + $(CC) $^ -o $@ $(LFLAGS) + @echo "done with $@" + +$(DESTBIN)/%: $(BIN)/%.$(ARCH) + @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi + rm -f $(DESTBIN)/$* + cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$* + +$(INSTALL): % : $(BIN)/%.$(ARCH) + +%.clean : + rm -f $(BIN)/$*.$(ARCH) + +%.install: + make $(DESTBIN)/$* + +# utilities ################################################# + +install: + for i in $(INSTALL); do make $$i.install; done + +clean: + rm -f $(BIN)/*.$(ARCH) + rm -f `find . -name "*.o"` + rm -f `find . -name "*~"` + rm -f `find . -name "#*"` Index: unk/archive/scripts/phase2-alt.c =================================================================== --- /trunk/archive/scripts/phase2-alt.c (revision 3702) +++ (revision ) @@ -1,307 +1,0 @@ -# include - -main (int argc, char **argv) { - - psArray *headers; - psMetadata *base, *site, *camera, *recipe; - - pmConfigLoadBasic (&site, &base); - - fits = psFitsAlloc (input); - headers = psFitsReadHeaderSet (fits); - - pmConfigLoadCamera (&camera, site, headers->data[0]); - pmConfigLoadRecipe (&recipe, camera, "PHASE2"); - - // check the validity of the given set of headers based on expectactions for this camera - pmCameraValidateHeaders (headers, camera); - - // construct a complete, empty (ie, no pixel data) FPA structure for this camera / header set - FPA = pmFPAfromHeader (headers, camera); - - /* ideally, the sequence of cells will be the same as the sequence of extensions */ - 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]; - - cellMD = psMetadataLookupPtr (cell->metadata, "CELL"); - Nextend = psMetadataLookupS32 (cellMD, "EXT.NUMBER"); - f = psFitsMoveExtNum(f, extNum, FALSE); // this should not move if we are already there - - for (k = 0; k < cell[0].readouts.n; k++) { - readout = cell[0].readouts[k]; - - // load the specific readout pixels. put this in a function? - cell[0].readouts[k] = pmReadoutLoad (cell[0].readouts[k], fits, cell, k); - bias_subtract_readout (readout, NULL, 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 = pmConfigLookupRegion (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); - -} - -// 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); -} - -bool pmConfigLoadSite (psMetadata **site, int *argc, char **argv) { - - int N; - char *sitefile; - - if ((N = psArgsGet (argc, argv, "-site"))) { - if (N == argc - 1) { - psError ("syntax error: -site (file)"); - return (FALSE); - } - psArgsRemove (&argc, argv, N); - sitefile = psStringCopy (argv[N]); - psArgsRemove (&argc, argv, N); - } - - if (sitefile == NULL) { - sitefile = psStringCopy ("site.cnf"); - } - *site = psMetadataParseConfig (NULL, sitefile, FALSE); - return (TRUE); -} - -bool pmConfigLoadCamera (psMetadata **camera, psMetadata *site, psMetadata *header) { - - int N; - char *cameraname = NULL; - char *camerafile = NULL; - - if (site == NULL) return (FALSE); - - // try to determine camera from command-line arguments - if ((N = psArgsGet (argc, argv, "-camera"))) { - if (N == argc - 1) { - psError ("syntax error: -camera (file)"); - return (FALSE); - } - psArgsRemove (&argc, argv, N); - cameraname = psStringCopy (argv[N]); - psArgsRemove (&argc, argv, N); - } - - // try to determine camera from header - if ((cameraname == NULL) && (header != NULL)) { - cameraname = pmCameraFromHeader (header); - } - - // try to determine camera from general configuration information - if (cameraname == NULL) { - cameraname = psMetadataLookupPtr (site, "CAMERA"); - } - - // if there is no source for camera info, fail - if (cameraname == NULL) { - return (FALSE); - } - - // find appropriate camera definition file - camerafile = psMetadataLookup (site, "CAMERA.%s", cameraname); - - // load the camera details for the specific camera - *camera = psMetadataParseConfig (NULL, camerafile, FALSE); - - return (TRUE); -} - -bool pmConfigLoadRecipe (psMetadata **recipe, psMetadata *camera, int *argc, char **argv, char *script) { - - int N; - char *recipename = NULL; - char *recipefile = NULL; - - if (camera == NULL) return (FALSE); - - // try to determine camera from command-line arguments - if ((N = psArgsGet (argc, argv, "-recipe"))) { - if (N == argc - 1) { - psError ("syntax error: -recipe (file)"); - return (FALSE); - } - psArgsRemove (&argc, argv, N); - recipefile = psStringCopy (argv[N]); - psArgsRemove (&argc, argv, N); - } - - // recipefile is defined in the camera metadata, if not specified on command line - if (recipefile == NULL) { - recipename = psAlloc (sizeof(script) + 8); - sprintf (recipename, "RECIPE.%s", script); - - // option 1: recipe is in the camera definition - recipefile = psMetadataLookup (camera, recipename); - psFree (recipename); - - // option 2: recipe is in the metadata database (equiv to a detrend file) - // this option requires a different API (header data, etc) - // recipefile = DetrendDatabaseLookup (recipename, site, camera, header); - } - - *recipe = psMetadataParseConfig (NULL, recipefile, FALSE); - psFree (recipefile); - - return (TRUE); -} Index: unk/archive/scripts/phase2.c =================================================================== --- /trunk/archive/scripts/phase2.c (revision 3702) +++ (revision ) @@ -1,364 +1,0 @@ -# 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]; - - fits = psFitsAlloc (input); - - // load only the headers from the given file - headers = psFitsReadHeaderSet (fits); - - // identify camera (may be specified on command line) - if (cameraname == NULL) { - cameraname = pmCameraFromHeader (psListGet (headers->list, 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) - - /* This will work, but it is not efficient: - we should work with the extensions in sequence */ - for (i = 0; i < FPA[0].chips.n; i++) { - chip = FPA[0].chips.data[i]; - - if (Base == CELL) { - f = psFitsMoveExtName(f, extname); - } - - for (j = 0; j < chip[0].cells.n; j++) { - cell = chip[0].cells.data[j]; - - if (Base == CELL) { - f = psFitsMoveExtName(f, extname); - } - - /* these need to work with the psFits pointer correctly - 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? - status = pmReadoutLoad (fits, cell, k); - - 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); - -} - - -bool pmConfigLoadBasic (psMetadata **site, psMetadata **base) { - - // 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); - -} - -bool pmConfigLoadCamera (psMetadata **camera, psMetadata *site, psMetadata *header) { - - if (site == NULL) return (FALSE); - - // identify camera (may be specified on command line) - cameraname = getargs ("-camera"); // pseudo-code, not a valid API - - // if there is no source for camera info, fail - if ((cameraname == NULL) && (header == NULL)) { - return (FALSE); - } - - // try to determine camera from header - if (cameraname == NULL) { - cameraname = pmCameraFromHeader (header); - } - - // find appropriate camera definition file - camerafile = psMetadataLookup (site, "CAMERA.%s", cameraname); - - // load the camera details for the specific camera - *camera = psMetadataParseConfig (NULL, camerafile, FALSE); - - return (TRUE); -} - -bool pmConfigLoadRecipe (psMetadata **recipe, psMetadata *camera, char *script) { - - char recipename[64]; - - if (site == NULL) return (FALSE); - if (camera == NULL) return (FALSE); - - // define command-line options - recipefile = getargs ("-recipe"); // pseudo-code, not a valid API - - // recipefile is defined in the camera metadata, if not specified on command line - if (recipefile == NULL) { - sprintf (recipename, "RECIPE.%s", script); - - // option 1: recipe is in the camera definition - recipefile = psMetadataLookup (camera, recipename); - - // option 2: recipe is in the metadata database (equiv to a detrend file) - // this option requires a different API (header data, etc) - recipefile = DetrendDatabaseLookup (recipename, site, camera, header); - } - - *recipe = psMetadataParseConfig (NULL, recipefile, FALSE); - return (TRUE); -} Index: unk/archive/scripts/phase2.pl =================================================================== --- /trunk/archive/scripts/phase2.pl (revision 3702) +++ (revision ) @@ -1,65 +1,0 @@ -#!/usr/bin/env perl - -sub usage { - # input is foo.fits - # outputs are outroot.im.fits (image) outroot.ob.fits (object list), outroot.log, etc - print STDERR "USAGE: phase2 (input) (outroot)\n"; - exit 2; -} - -if (@ARGV != 4) { &usage; } - -# base file needs to be hardwired, from environment, or based on program name -$base = psMetadataParseConfig (NULL, basefile, FALSE); - -# probably defined in the base file? -$site = psMetadataParseConfig (NULL, sitefile, FALSE); - -# this should be optionally specified as part of the command-line options -# alternatively, this should be derived from the input file header data -$camera = psMetadataParseConfig (NULL, camerafile, FALSE); - -# we need a way to define an optional alternative recipe -$recipe = psMetadataParseConfig (NULL, recipefile, FALSE); # defined in the camera file - -# load only the headers from the given file input a list of headers -@headers = psFITSReadHeaderSet ($input); - -# check the validity of the given set of headers based on expectactions for this camera -# ** where does camera come from? -# ** how do we do this if the camera is derived from the header info? -$status = pmCameraValidateHeaders ($headers, $camera); - -# this gets a bit tenuous: what do we do with the collection? - -# 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); - -### now we work on individual readouts: - -# we have the FPA, already validated -# as part of the FPA, we have access to a specified READOUT and we know the corresponding EXTNAME - -# load in the complete readout: -$image = psImageReadSection (NULL, 0, 0, 0, 0, 0, $extname, $Nreadout, $input); - -# load in the complete corresponding bias image (this should be done once per cell) -# in this form, we are implying that the bias image is guaranteed to have the corresponding EXTNAME -$bias = psImageReadSection (NULL, 0, 0, 0, 0, 0, $extname, 0, $biasname); - -# overscan options: -# - overscan is defined by the recipe (BIAS.OVERSCAN has form [Xs,Xe:Ys:Ye]) -# - overscan is defined by the header (BIAS.OVERSCAN has value BIASSEC) -# - overscan is not used (BIAS.OVERSCAN has value NONE) - -# bias section from header: -$ = psMetadataLookup ($readout->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: unk/archive/scripts/phase4.c =================================================================== --- /trunk/archive/scripts/phase4.c (revision 3702) +++ (revision ) @@ -1,212 +1,0 @@ -#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: unk/archive/scripts/pmFPAio.c =================================================================== --- /trunk/archive/scripts/pmFPAio.c (revision 3702) +++ (revision ) @@ -1,20 +1,0 @@ -# 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? Index: /trunk/archive/scripts/scripts/phase2.pl =================================================================== --- /trunk/archive/scripts/scripts/phase2.pl (revision 3703) +++ /trunk/archive/scripts/scripts/phase2.pl (revision 3703) @@ -0,0 +1,65 @@ +#!/usr/bin/env perl + +sub usage { + # input is foo.fits + # outputs are outroot.im.fits (image) outroot.ob.fits (object list), outroot.log, etc + print STDERR "USAGE: phase2 (input) (outroot)\n"; + exit 2; +} + +if (@ARGV != 4) { &usage; } + +# base file needs to be hardwired, from environment, or based on program name +$base = psMetadataParseConfig (NULL, basefile, FALSE); + +# probably defined in the base file? +$site = psMetadataParseConfig (NULL, sitefile, FALSE); + +# this should be optionally specified as part of the command-line options +# alternatively, this should be derived from the input file header data +$camera = psMetadataParseConfig (NULL, camerafile, FALSE); + +# we need a way to define an optional alternative recipe +$recipe = psMetadataParseConfig (NULL, recipefile, FALSE); # defined in the camera file + +# load only the headers from the given file input a list of headers +@headers = psFITSReadHeaderSet ($input); + +# check the validity of the given set of headers based on expectactions for this camera +# ** where does camera come from? +# ** how do we do this if the camera is derived from the header info? +$status = pmCameraValidateHeaders ($headers, $camera); + +# this gets a bit tenuous: what do we do with the collection? + +# 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); + +### now we work on individual readouts: + +# we have the FPA, already validated +# as part of the FPA, we have access to a specified READOUT and we know the corresponding EXTNAME + +# load in the complete readout: +$image = psImageReadSection (NULL, 0, 0, 0, 0, 0, $extname, $Nreadout, $input); + +# load in the complete corresponding bias image (this should be done once per cell) +# in this form, we are implying that the bias image is guaranteed to have the corresponding EXTNAME +$bias = psImageReadSection (NULL, 0, 0, 0, 0, 0, $extname, 0, $biasname); + +# overscan options: +# - overscan is defined by the recipe (BIAS.OVERSCAN has form [Xs,Xe:Ys:Ye]) +# - overscan is defined by the header (BIAS.OVERSCAN has value BIASSEC) +# - overscan is not used (BIAS.OVERSCAN has value NONE) + +# bias section from header: +$ = psMetadataLookup ($readout->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/src/phase2-alt.c =================================================================== --- /trunk/archive/scripts/src/phase2-alt.c (revision 3703) +++ /trunk/archive/scripts/src/phase2-alt.c (revision 3703) @@ -0,0 +1,307 @@ +# include + +main (int argc, char **argv) { + + psArray *headers; + psMetadata *base, *site, *camera, *recipe; + + pmConfigLoadBasic (&site, &base); + + fits = psFitsAlloc (input); + headers = psFitsReadHeaderSet (fits); + + pmConfigLoadCamera (&camera, site, headers->data[0]); + pmConfigLoadRecipe (&recipe, camera, "PHASE2"); + + // check the validity of the given set of headers based on expectactions for this camera + pmCameraValidateHeaders (headers, camera); + + // construct a complete, empty (ie, no pixel data) FPA structure for this camera / header set + FPA = pmFPAfromHeader (headers, camera); + + /* ideally, the sequence of cells will be the same as the sequence of extensions */ + 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]; + + cellMD = psMetadataLookupPtr (cell->metadata, "CELL"); + Nextend = psMetadataLookupS32 (cellMD, "EXT.NUMBER"); + f = psFitsMoveExtNum(f, extNum, FALSE); // this should not move if we are already there + + for (k = 0; k < cell[0].readouts.n; k++) { + readout = cell[0].readouts[k]; + + // load the specific readout pixels. put this in a function? + cell[0].readouts[k] = pmReadoutLoad (cell[0].readouts[k], fits, cell, k); + bias_subtract_readout (readout, NULL, 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 = pmConfigLookupRegion (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); + +} + +// 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); +} + +bool pmConfigLoadSite (psMetadata **site, int *argc, char **argv) { + + int N; + char *sitefile; + + if ((N = psArgsGet (argc, argv, "-site"))) { + if (N == argc - 1) { + psError ("syntax error: -site (file)"); + return (FALSE); + } + psArgsRemove (&argc, argv, N); + sitefile = psStringCopy (argv[N]); + psArgsRemove (&argc, argv, N); + } + + if (sitefile == NULL) { + sitefile = psStringCopy ("site.cnf"); + } + *site = psMetadataParseConfig (NULL, sitefile, FALSE); + return (TRUE); +} + +bool pmConfigLoadCamera (psMetadata **camera, psMetadata *site, psMetadata *header) { + + int N; + char *cameraname = NULL; + char *camerafile = NULL; + + if (site == NULL) return (FALSE); + + // try to determine camera from command-line arguments + if ((N = psArgsGet (argc, argv, "-camera"))) { + if (N == argc - 1) { + psError ("syntax error: -camera (file)"); + return (FALSE); + } + psArgsRemove (&argc, argv, N); + cameraname = psStringCopy (argv[N]); + psArgsRemove (&argc, argv, N); + } + + // try to determine camera from header + if ((cameraname == NULL) && (header != NULL)) { + cameraname = pmCameraFromHeader (header); + } + + // try to determine camera from general configuration information + if (cameraname == NULL) { + cameraname = psMetadataLookupPtr (site, "CAMERA"); + } + + // if there is no source for camera info, fail + if (cameraname == NULL) { + return (FALSE); + } + + // find appropriate camera definition file + camerafile = psMetadataLookup (site, "CAMERA.%s", cameraname); + + // load the camera details for the specific camera + *camera = psMetadataParseConfig (NULL, camerafile, FALSE); + + return (TRUE); +} + +bool pmConfigLoadRecipe (psMetadata **recipe, psMetadata *camera, int *argc, char **argv, char *script) { + + int N; + char *recipename = NULL; + char *recipefile = NULL; + + if (camera == NULL) return (FALSE); + + // try to determine camera from command-line arguments + if ((N = psArgsGet (argc, argv, "-recipe"))) { + if (N == argc - 1) { + psError ("syntax error: -recipe (file)"); + return (FALSE); + } + psArgsRemove (&argc, argv, N); + recipefile = psStringCopy (argv[N]); + psArgsRemove (&argc, argv, N); + } + + // recipefile is defined in the camera metadata, if not specified on command line + if (recipefile == NULL) { + recipename = psAlloc (sizeof(script) + 8); + sprintf (recipename, "RECIPE.%s", script); + + // option 1: recipe is in the camera definition + recipefile = psMetadataLookup (camera, recipename); + psFree (recipename); + + // option 2: recipe is in the metadata database (equiv to a detrend file) + // this option requires a different API (header data, etc) + // recipefile = DetrendDatabaseLookup (recipename, site, camera, header); + } + + *recipe = psMetadataParseConfig (NULL, recipefile, FALSE); + psFree (recipefile); + + return (TRUE); +} Index: /trunk/archive/scripts/src/phase2.c =================================================================== --- /trunk/archive/scripts/src/phase2.c (revision 3703) +++ /trunk/archive/scripts/src/phase2.c (revision 3703) @@ -0,0 +1,364 @@ +# 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]; + + fits = psFitsAlloc (input); + + // load only the headers from the given file + headers = psFitsReadHeaderSet (fits); + + // identify camera (may be specified on command line) + if (cameraname == NULL) { + cameraname = pmCameraFromHeader (psListGet (headers->list, 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) + + /* This will work, but it is not efficient: + we should work with the extensions in sequence */ + for (i = 0; i < FPA[0].chips.n; i++) { + chip = FPA[0].chips.data[i]; + + if (Base == CELL) { + f = psFitsMoveExtName(f, extname); + } + + for (j = 0; j < chip[0].cells.n; j++) { + cell = chip[0].cells.data[j]; + + if (Base == CELL) { + f = psFitsMoveExtName(f, extname); + } + + /* these need to work with the psFits pointer correctly + 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? + status = pmReadoutLoad (fits, cell, k); + + 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); + +} + + +bool pmConfigLoadBasic (psMetadata **site, psMetadata **base) { + + // 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); + +} + +bool pmConfigLoadCamera (psMetadata **camera, psMetadata *site, psMetadata *header) { + + if (site == NULL) return (FALSE); + + // identify camera (may be specified on command line) + cameraname = getargs ("-camera"); // pseudo-code, not a valid API + + // if there is no source for camera info, fail + if ((cameraname == NULL) && (header == NULL)) { + return (FALSE); + } + + // try to determine camera from header + if (cameraname == NULL) { + cameraname = pmCameraFromHeader (header); + } + + // find appropriate camera definition file + camerafile = psMetadataLookup (site, "CAMERA.%s", cameraname); + + // load the camera details for the specific camera + *camera = psMetadataParseConfig (NULL, camerafile, FALSE); + + return (TRUE); +} + +bool pmConfigLoadRecipe (psMetadata **recipe, psMetadata *camera, char *script) { + + char recipename[64]; + + if (site == NULL) return (FALSE); + if (camera == NULL) return (FALSE); + + // define command-line options + recipefile = getargs ("-recipe"); // pseudo-code, not a valid API + + // recipefile is defined in the camera metadata, if not specified on command line + if (recipefile == NULL) { + sprintf (recipename, "RECIPE.%s", script); + + // option 1: recipe is in the camera definition + recipefile = psMetadataLookup (camera, recipename); + + // option 2: recipe is in the metadata database (equiv to a detrend file) + // this option requires a different API (header data, etc) + recipefile = DetrendDatabaseLookup (recipename, site, camera, header); + } + + *recipe = psMetadataParseConfig (NULL, recipefile, FALSE); + return (TRUE); +} Index: /trunk/archive/scripts/src/phase4.c =================================================================== --- /trunk/archive/scripts/src/phase4.c (revision 3703) +++ /trunk/archive/scripts/src/phase4.c (revision 3703) @@ -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/src/pmFPAio.c =================================================================== --- /trunk/archive/scripts/src/pmFPAio.c (revision 3703) +++ /trunk/archive/scripts/src/pmFPAio.c (revision 3703) @@ -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?