Index: /trunk/psModules/src/config/pmConfig.c =================================================================== --- /trunk/psModules/src/config/pmConfig.c (revision 14504) +++ /trunk/psModules/src/config/pmConfig.c (revision 14505) @@ -4,6 +4,6 @@ * @author EAM (IfA) * - * @version $Revision: 1.100 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-07-27 00:11:03 $ + * @version $Revision: 1.101 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-15 20:21:18 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -552,5 +552,9 @@ if (newFile) { cameraFile = newFile; - } + } else { + psLogMsg ("psModules.config", PS_LOG_INFO, + "%s is not listed in the site CAMERAS list,\n" + "trying as a camera configuration file\n", cameraFile); + } } else { psError(PS_ERR_IO, false, "Unable to find CAMERAS in site configuration.\n"); Index: /trunk/psModules/src/detrend/Makefile.am =================================================================== --- /trunk/psModules/src/detrend/Makefile.am (revision 14504) +++ /trunk/psModules/src/detrend/Makefile.am (revision 14505) @@ -10,4 +10,5 @@ pmNonLinear.c \ pmBias.c \ + pmOverscan.c \ pmDetrendDB.c \ pmShutterCorrection.c \ @@ -23,4 +24,5 @@ pmNonLinear.h \ pmBias.h \ + pmOverscan.h \ pmDetrendDB.h \ pmShutterCorrection.h \ Index: /trunk/psModules/src/detrend/pmBias.c =================================================================== --- /trunk/psModules/src/detrend/pmBias.c (revision 14504) +++ /trunk/psModules/src/detrend/pmBias.c (revision 14505) @@ -15,43 +15,10 @@ #include "pmFPACalibration.h" +#include "pmOverscan.h" #include "pmBias.h" -#define SMOOTH_NSIGMA 4.0 // Number of Gaussian sigma the smoothing kernel extends - - -static void overscanOptionsFree(pmOverscanOptions *options) -{ - psFree(options->stat); - psFree(options->poly); - psFree(options->spline); -} - -pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat, - int boxcar, float gauss) -{ - pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); - psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree); - - // Inputs - opts->single = single; - opts->fitType = fitType; - opts->order = order; - opts->stat = psMemIncrRefCounter(stat); - - // Smoothing - opts->boxcar = boxcar; - opts->gauss = gauss; - - // Outputs - opts->poly = NULL; - opts->spline = NULL; - - return opts; -} - - -// subtractFrame(): this routine will take as input a readout for the input image and a readout for the bias +// pmBiasSubtractFrame(): this routine will take as input a readout for the input image and a readout for the bias // image. The bias image is subtracted in place from the input image. -static bool subtractFrame(pmReadout *in, // Input readout +bool pmBiasSubtractFrame(pmReadout *in, // Input readout const pmReadout *sub, // Readout to be subtracted from input float scale // Scale to apply before subtracting @@ -125,345 +92,32 @@ } -// Produce an overscan vector from an array of pixels -static psVector *overscanVector(float *chi2, // chi^2 from fit - pmOverscanOptions *overscanOpts, // Overscan options - const psArray *pixels, // Array of vectors containing the pixel values - psStats *myStats // Statistic to use in reducing the overscan - ) -{ - assert(overscanOpts); - assert(pixels); - assert(myStats); - - psStatsOptions statistic = psStatsSingleOption(myStats->options); // Statistic to use - assert(statistic != 0); - - // Reduce the overscans - psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row - psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate - psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting - - for (int i = 0; i < pixels->n; i++) { - psVector *values = pixels->data[i]; // Vector with overscan values - if (values->n > 0) { - mask->data.U8[i] = 0; - ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1] - psVectorStats(myStats, values, NULL, NULL, 0); - reduced->data.F32[i] = psStatsGetValue(myStats, statistic); - } else if (overscanOpts->fitType == PM_FIT_NONE) { - psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the " - "image, and no fit is requested.\n"); - return NULL; - } else { - // We'll fit this one out - mask->data.U8[i] = 1; - } - } - - // Smooth the reduced vector - if (overscanOpts->boxcar > 0) { - psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar); // Smoothed vector - psFree(reduced); - reduced = smoothed; - } - if (isfinite(overscanOpts->gauss) && overscanOpts->gauss > 0) { - if (overscanOpts->boxcar > 0) { - psWarning("Gaussian smoothing the boxcar smoothed overscan --- you asked for it."); - } - psVector *smoothed = psVectorSmooth(NULL, reduced, overscanOpts->gauss, SMOOTH_NSIGMA); - psFree(reduced); - reduced = smoothed; - } - - // Fit the overscan, if required - psVector *fitted; // Fitted overscan values - switch (overscanOpts->fitType) { - case PM_FIT_NONE: - // No fitting --- that's easy. - fitted = psMemIncrRefCounter(reduced); - break; - case PM_FIT_POLY_ORD: - if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || - overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) { - psFree(overscanOpts->poly); - overscanOpts->poly = NULL; - } - if (! overscanOpts->poly) { - overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order); - } - psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); - fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); - break; - case PM_FIT_POLY_CHEBY: - if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || - overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) { - psFree(overscanOpts->poly); - overscanOpts->poly = NULL; - } - if (! overscanOpts->poly) { - overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order); - } - psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); - fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); - break; - case PM_FIT_SPLINE: - // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an - // input spline - overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); - fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate); - break; - default: - psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); - return NULL; - break; - } - - if (chi2) { - *chi2 = 0.0; // chi^2 (sort of) - for (int i = 0; i < reduced->n; i++) { - *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]); - } - } - - psFree(reduced); - psFree(ordinate); - psFree(mask); - - return fitted; -} - -bool pmBiasSubtract(pmReadout *inRO, pmOverscanOptions *overscanOpts, - const pmReadout *biasRO, const pmReadout *darkRO, const pmFPAview *view) +bool pmBiasSubtract(pmReadout *in, pmOverscanOptions *overscanOpts, + const pmReadout *bias, const pmReadout *dark, const pmFPAview *view) { psTrace("psModules.detrend", 4, "---- pmBiasSubtract() begin ----\n"); - PS_ASSERT_PTR_NON_NULL(inRO, false); - PS_ASSERT_IMAGE_NON_NULL(inRO->image, false); - PS_ASSERT_IMAGE_TYPE(inRO->image, PS_TYPE_F32, false); - if (biasRO) { - PS_ASSERT_IMAGE_NON_NULL(biasRO->image, false); - PS_ASSERT_IMAGE_TYPE(biasRO->image, PS_TYPE_F32, false); + PS_ASSERT_PTR_NON_NULL(in, false); + PS_ASSERT_IMAGE_NON_NULL(in->image, false); + PS_ASSERT_IMAGE_TYPE(in->image, PS_TYPE_F32, false); + if (bias) { + PS_ASSERT_IMAGE_NON_NULL(bias->image, false); + PS_ASSERT_IMAGE_TYPE(bias->image, PS_TYPE_F32, false); } - if (darkRO) { + if (dark) { PS_ASSERT_PTR_NON_NULL(view, false); - PS_ASSERT_IMAGE_NON_NULL(darkRO->image, false); - PS_ASSERT_IMAGE_TYPE(darkRO->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(dark->image, false); + PS_ASSERT_IMAGE_TYPE(dark->image, PS_TYPE_F32, false); } - pmHDU *hdu = pmHDUFromReadout(inRO); // HDU of interest - psImage *image = inRO->image; // The input image + pmHDU *hdu = pmHDUFromReadout(in); // HDU of interest - // Overscan processing - if (overscanOpts) { - // Check for an unallowable pmFit. - if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD && - overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) { - psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d). Returning original image.\n", - overscanOpts->fitType); - return false; - } - - psList *overscans = inRO->bias; // List of the overscan images - - psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use - if (statistic == 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", - overscanOpts->stat); - return false; - } - psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original - - psString comment = NULL; // Comment to add - psStringAppend(&comment, "Subtracting overscan (stat %x; type %x; order %d)", - statistic, overscanOpts->fitType, overscanOpts->order); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, - comment, ""); - psFree(comment); - - // Reduce all overscan pixels to a single value - if (overscanOpts->single) { - psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - int index = pixels->n; // Index - pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols); - pixels->n += overscan->numRows * overscan->numCols; - for (int i = 0; i < overscan->numRows; i++) { - memcpy(&pixels->data.F32[index], overscan->data.F32[i], - overscan->numCols * sizeof(psF32)); - index += overscan->numCols; - } - } - psFree(iter); - - (void)psVectorStats(stats, pixels, NULL, NULL, 0); - psFree(pixels); - double reduced = psStatsGetValue(stats, statistic); // Result of statistics - - psString comment = NULL; // Comment to add - psStringAppend(&comment, "Overscan value: %f", reduced); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", - reduced); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); - - (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32)); - } else { - // We do the regular overscan subtraction - bool readRows = psMetadataLookupBool(NULL, inRO->parent->concepts, - "CELL.READDIR"); // Read direction - float chi2 = NAN; // chi^2 from fit - - if (readRows) { - // The read direction is rows - psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels - for (int i = 0; i < pixels->n; i++) { - pixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); - } - - // Pull the pixels out into the vectors - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - // the overscan and image might not be aligned. pixels->data represents - // the image row pixels. - int diff = overscan->row0 - image->row0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) { - int j = i - diff; - // i is row on image - // j is row on overscan - psVector *values = pixels->data[i]; - int index = values->n; // Index in the vector - values = psVectorRealloc(values, values->n + overscan->numCols); - values->n += overscan->numCols; - memcpy(&values->data.F32[index], overscan->data.F32[j], - overscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); - index += overscan->numCols; - pixels->data[i] = values; // Update the pointer in case it's moved - } - } - psFree(iter); - - // Reduce the overscans - psVector *reduced = overscanVector(&chi2, overscanOpts, pixels, stats); - psFree(pixels); - if (! reduced) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); - return false; - } - - // Subtract row by row - for (int i = 0; i < image->numRows; i++) { - for (int j = 0; j < image->numCols; j++) { - image->data.F32[i][j] -= reduced->data.F32[i]; - } - } - psFree(reduced); - - } else { - // The read direction is columns - psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels - for (int i = 0; i < pixels->n; i++) { - psVector *values = psVectorAlloc(0, PS_TYPE_F32); - pixels->data[i] = values; - } - - // Pull the pixels out into the vectors - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - // the overscan and image might not be aligned. pixels->data represents - // the image row pixels. - int diff = overscan->col0 - image->col0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) { - int iFixed = i - diff; - // i is column on image - // iFixed is column on overscan - psVector *values = pixels->data[i]; - int index = values->n; // Index in the vector - values = psVectorRealloc(values, values->n + overscan->numRows); - for (int j = 0; j < overscan->numRows; j++) { - values->data.F32[index++] = overscan->data.F32[j][iFixed]; - } - values->n += overscan->numRows; - pixels->data[i] = values; // Update the pointer in case it's moved - } - } - psFree(iter); - - // Reduce the overscans - psVector *reduced = overscanVector(&chi2, overscanOpts, pixels, stats); - psFree(pixels); - if (! reduced) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); - return false; - } - - // Subtract column by column - for (int i = 0; i < image->numCols; i++) { - for (int j = 0; j < image->numRows; j++) { - image->data.F32[j][i] -= reduced->data.F32[i]; - } - } - psFree(reduced); - } - - switch (overscanOpts->fitType) { - case PM_FIT_POLY_ORD: - case PM_FIT_POLY_CHEBY: { - psString comment = NULL; // Comment to add - psStringAppend(&comment, "Overscan fit (chi2: %.2f): ", chi2); - psPolynomial1D *poly = overscanOpts->poly; // The polynomial - for (int i = 0; i < poly->nX; i++) { - psStringAppend(&comment, "%.1f ", poly->coeff[i]); - } - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, - "Overscan value", poly->coeff[0]); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, - "Overscan stdev", poly->coeffErr[0]); - break; - } - case PM_FIT_SPLINE: { - psSpline1D *spline = overscanOpts->spline; // The spline - for (int i = 0; i < spline->n; i++) { - psStringAppend(&comment, "Overscan fit (chi2: %.2f) %d:", chi2, i); - psPolynomial1D *poly = spline->spline[i]; // i-th polynomial - for (int j = 0; j < poly->nX; j++) { - psStringAppend(&comment, "%.1f ", poly->coeff[i]); - } - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, - comment, ""); - psFree(comment); - } - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, - "Overscan value", NAN); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, - "Overscan stdev", NAN); - break; - } - case PM_FIT_NONE: - break; - default: - psAbort("Should never get here!!!\n"); - } - - - } - psFree(stats); - } // End of overscan subtraction + if (!pmOverscanSubtract (in, overscanOpts)) { + return false; + } // Bias frame subtraction - if (biasRO) { - psVector *md5 = psImageMD5(biasRO->image); // md5 hash + if (bias) { + psVector *md5 = psImageMD5(bias->image); // md5 hash psString md5string = psMD5toString(md5); // String psFree(md5); @@ -473,13 +127,13 @@ psFree(md5string); - if (!subtractFrame(inRO, biasRO, 1.0)) { + if (!pmBiasSubtractFrame(in, bias, 1.0)) { return false; } } - if (darkRO) { + if (dark) { // Get the scaling - float inTime = psMetadataLookupF32(NULL, inRO->parent->concepts, "CELL.DARKTIME"); - float darkTime = psMetadataLookupF32(NULL, darkRO->parent->concepts, "CELL.DARKTIME"); + float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME"); + float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME"); if (isnan(inTime) || isnan(darkTime)) { psError(PS_ERR_UNKNOWN, false, "Unable to determine dark scaling."); @@ -488,5 +142,5 @@ float darkNorm = 1.0; - float inNorm = pmFPADarkNorm(inRO->parent->parent->parent, view, inTime); + float inNorm = pmFPADarkNorm(in->parent->parent->parent, view, inTime); // if we have a normalized dark exposure, we simply multiply the master by inNorm. if @@ -495,5 +149,5 @@ if (darkTime != 1.0) { - darkNorm = pmFPADarkNorm(darkRO->parent->parent->parent, view, darkTime); + darkNorm = pmFPADarkNorm(dark->parent->parent->parent, view, darkTime); } @@ -505,5 +159,5 @@ float scale = inNorm / darkNorm;// Scaling to apply to dark exposure - psVector *md5 = psImageMD5(darkRO->image); // md5 hash + psVector *md5 = psImageMD5(dark->image); // md5 hash psString md5string = psMD5toString(md5); // String psFree(md5); @@ -513,5 +167,5 @@ psFree(md5string); - if (!subtractFrame(inRO, darkRO, scale)) { + if (!pmBiasSubtractFrame(in, dark, scale)) { return false; } Index: /trunk/psModules/src/detrend/pmBias.h =================================================================== --- /trunk/psModules/src/detrend/pmBias.h (revision 14504) +++ /trunk/psModules/src/detrend/pmBias.h (revision 14505) @@ -5,53 +5,14 @@ * @author Paul Price, IfA * - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-06-20 20:45:00 $ + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-15 20:21:18 $ * Copyright 2004--2006 Institute for Astronomy, University of Hawaii */ -#ifndef PM__BIAS_H -#define PM__BIAS_H +#ifndef PM_BIAS_H +#define PM_BIAS_H /// @addtogroup detrend Detrend Creation and Application /// @{ - -/// Type of fit to perform -typedef enum { - PM_FIT_NONE, ///< No fit - PM_FIT_POLY_ORD, ///< Fit ordinary polynomial - PM_FIT_POLY_CHEBY, ///< Fit Chebyshev polynomial - PM_FIT_SPLINE ///< Fit cubic splines -} pmFit; - -/// Options for overscan subtraction -/// -/// The overscan subtraction may be performed by reducing all overscan regions to a single value (e.g., if -/// there is no structure); or the overscan may be fit perpendicular to the read direction (usually the -/// columns) with a particular functional form; or a single value may be subtracted for each read/scan without -/// fitting (if the structure defies characterisation). In any case, statistics are required to reduce -/// multiple values to a single value (either for the scan, or for the entire overscan regions). -typedef struct -{ - // Inputs - bool single; ///< Reduce all overscan regions to a single value? - pmFit fitType; ///< Type of fit to overscan - unsigned int order; ///< Order of polynomial, or number of spline pieces - psStats *stat; ///< Statistic to use when reducing the minor direction - int boxcar; ///< Boxcar smoothing radius - float gauss; ///< Gaussian smoothing sigma - // Outputs - psPolynomial1D *poly; ///< Result of polynomial fit - psSpline1D *spline; ///< Result of spline fit -} -pmOverscanOptions; - -/// Allocator for overscan options -pmOverscanOptions *pmOverscanOptionsAlloc(bool single, ///< Reduce all overscan regions to a single value? - pmFit fitType, ///< Type of fit to overscan - unsigned int order, ///< Order of polynomial, or number of splines - psStats *stat, ///< Statistic to use - int boxcar, ///< Boxcar smoothing radius - float gauss ///< Gaussian smoothing sigma - ); /// Subtract the overscan, bias and/or dark @@ -66,4 +27,11 @@ ); +// pmBiasSubtractFrame(): this routine will take as input a readout for the input image and a readout for the bias +// image. The bias image is subtracted in place from the input image. +bool pmBiasSubtractFrame(pmReadout *in, // Input readout + const pmReadout *sub, // Readout to be subtracted from input + float scale // Scale to apply before subtracting + ); + /// @} #endif Index: /trunk/psModules/src/detrend/pmOverscan.c =================================================================== --- /trunk/psModules/src/detrend/pmOverscan.c (revision 14505) +++ /trunk/psModules/src/detrend/pmOverscan.c (revision 14505) @@ -0,0 +1,398 @@ +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmHDUUtils.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPACalibration.h" + +#include "pmOverscan.h" + +#define SMOOTH_NSIGMA 4.0 // Number of Gaussian sigma the smoothing kernel extends + +static void pmOverscanOptionsFree(pmOverscanOptions *options) +{ + psFree(options->stat); + psFree(options->poly); + psFree(options->spline); +} + +pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat, + int boxcar, float gauss) +{ + pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); + psMemSetDeallocator(opts, (psFreeFunc)pmOverscanOptionsFree); + + // Inputs + opts->single = single; + opts->constant = false; + opts->fitType = fitType; + opts->order = order; + opts->stat = psMemIncrRefCounter(stat); + + // Smoothing + opts->boxcar = boxcar; + opts->gauss = gauss; + + // Outputs + opts->poly = NULL; + opts->spline = NULL; + + return opts; +} + +// Produce an overscan vector from an array of pixels +psVector *pmOverscanVector(float *chi2, // chi^2 from fit + pmOverscanOptions *overscanOpts, // Overscan options + const psArray *pixels, // Array of vectors containing the pixel values + psStats *myStats // Statistic to use in reducing the overscan + ) +{ + assert(overscanOpts); + assert(pixels); + assert(myStats); + + psStatsOptions statistic = psStatsSingleOption(myStats->options); // Statistic to use + assert(statistic != 0); + + // Reduce the overscans + psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row + psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate + psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting + + for (int i = 0; i < pixels->n; i++) { + psVector *values = pixels->data[i]; // Vector with overscan values + if (values->n > 0) { + mask->data.U8[i] = 0; + ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1] + psVectorStats(myStats, values, NULL, NULL, 0); + reduced->data.F32[i] = psStatsGetValue(myStats, statistic); + } else if (overscanOpts->fitType == PM_FIT_NONE) { + psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the " + "image, and no fit is requested.\n"); + return NULL; + } else { + // We'll fit this one out + mask->data.U8[i] = 1; + } + } + + // Smooth the reduced vector + if (overscanOpts->boxcar > 0) { + psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar); // Smoothed vector + psFree(reduced); + reduced = smoothed; + } + if (isfinite(overscanOpts->gauss) && overscanOpts->gauss > 0) { + if (overscanOpts->boxcar > 0) { + psWarning("Gaussian smoothing the boxcar smoothed overscan --- you asked for it."); + } + psVector *smoothed = psVectorSmooth(NULL, reduced, overscanOpts->gauss, SMOOTH_NSIGMA); + psFree(reduced); + reduced = smoothed; + } + + // Fit the overscan, if required + psVector *fitted; // Fitted overscan values + switch (overscanOpts->fitType) { + case PM_FIT_NONE: + // No fitting --- that's easy. + fitted = psMemIncrRefCounter(reduced); + break; + case PM_FIT_POLY_ORD: + if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || + overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) { + psFree(overscanOpts->poly); + overscanOpts->poly = NULL; + } + if (! overscanOpts->poly) { + overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order); + } + psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); + fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); + break; + case PM_FIT_POLY_CHEBY: + if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || + overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) { + psFree(overscanOpts->poly); + overscanOpts->poly = NULL; + } + if (! overscanOpts->poly) { + overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order); + } + psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); + fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); + break; + case PM_FIT_SPLINE: + // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an + // input spline + overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); + fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); + return NULL; + break; + } + + if (chi2) { + *chi2 = 0.0; // chi^2 (sort of) + for (int i = 0; i < reduced->n; i++) { + *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]); + } + } + + psFree(reduced); + psFree(ordinate); + psFree(mask); + + return fitted; +} + +bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2) { + + psString comment = NULL; // Comment to add + + switch (overscanOpts->fitType) { + case PM_FIT_POLY_ORD: + case PM_FIT_POLY_CHEBY: { + psStringAppend(&comment, "Overscan fit (chi2: %.2f): ", chi2); + psPolynomial1D *poly = overscanOpts->poly; // The polynomial + for (int i = 0; i < poly->nX; i++) { + psStringAppend(&comment, "%.1f ", poly->coeff[i]); + } + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + comment = NULL; + + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, + "Overscan value", poly->coeff[0]); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, + "Overscan stdev", poly->coeffErr[0]); + break; + } + case PM_FIT_SPLINE: { + psSpline1D *spline = overscanOpts->spline; // The spline + for (int i = 0; i < spline->n; i++) { + psStringAppend(&comment, "Overscan fit (chi2: %.2f) %d:", chi2, i); + psPolynomial1D *poly = spline->spline[i]; // i-th polynomial + for (int j = 0; j < poly->nX; j++) { + psStringAppend(&comment, "%.1f ", poly->coeff[i]); + } + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, + comment, ""); + psFree(comment); + comment = NULL; + } + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, + "Overscan value", NAN); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, + "Overscan stdev", NAN); + break; + } + case PM_FIT_NONE: + break; + default: + psAbort("Should never get here!!!\n"); + } + return true; +} + +bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts) { + + assert (input); + + if (overscanOpts == NULL) return true; // no overscan subtraction requested + + pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest + psImage *image = input->image; + + // check for 'soft bias' (simple, fixed offset to be subtracted) + if (overscanOpts->constant) { + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", + overscanOpts->value); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); + + (void)psBinaryOp(input->image, input->image, "-", psScalarAlloc((float)overscanOpts->value, PS_TYPE_F32)); + + return true; + } + + // we are performing a statitical analysis of the overscan region + + // Check for an unallowable pmFit. + if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD && + overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) { + psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d). Returning original image.\n", + overscanOpts->fitType); + return false; + } + + psList *overscans = input->bias; // List of the overscan images + + psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use + if (statistic == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", + overscanOpts->stat); + return false; + } + psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original + + psString comment = NULL; // Comment to add + psStringAppend(&comment, "Subtracting overscan (stat %x; type %x; order %d)", + statistic, overscanOpts->fitType, overscanOpts->order); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, + comment, ""); + psFree(comment); + + // Reduce all overscan pixels to a single value + if (overscanOpts->single) { + psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) { + int index = pixels->n; // Index + pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols); + pixels->n += overscan->numRows * overscan->numCols; + for (int i = 0; i < overscan->numRows; i++) { + memcpy(&pixels->data.F32[index], overscan->data.F32[i], + overscan->numCols * sizeof(psF32)); + index += overscan->numCols; + } + } + psFree(iter); + + (void)psVectorStats(stats, pixels, NULL, NULL, 0); + psFree(pixels); + double reduced = psStatsGetValue(stats, statistic); // Result of statistics + + psString comment = NULL; // Comment to add + psStringAppend(&comment, "Overscan value: %f", reduced); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", + reduced); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); + + (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32)); + return true; + } + + // We are performing a row-by-row overscan subtraction + bool readRows = psMetadataLookupBool(NULL, input->parent->concepts, "CELL.READDIR"); // Read direction + float chi2 = NAN; // chi^2 from fit + + // adjust operation depending on the read direction : need to re-org pixels for columns + if (readRows) { + // The read direction is rows + psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels + for (int i = 0; i < pixels->n; i++) { + pixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); + } + + // Pull the pixels out into the vectors + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) { + // the overscan and image might not be aligned. pixels->data represents + // the image row pixels. + int diff = overscan->row0 - image->row0; // Offset between the two regions + for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) { + int j = i - diff; + // i is row on image + // j is row on overscan + psVector *values = pixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + overscan->numCols); + values->n += overscan->numCols; + memcpy(&values->data.F32[index], overscan->data.F32[j], + overscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + index += overscan->numCols; + pixels->data[i] = values; // Update the pointer in case it's moved + } + } + psFree(iter); + + // Reduce the overscans + psVector *reduced = pmOverscanVector(&chi2, overscanOpts, pixels, stats); + psFree(pixels); + if (! reduced) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); + return false; + } + + // Subtract row by row + for (int i = 0; i < image->numRows; i++) { + for (int j = 0; j < image->numCols; j++) { + image->data.F32[i][j] -= reduced->data.F32[i]; + } + } + psFree(reduced); + } else { + // The read direction is columns + psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels + for (int i = 0; i < pixels->n; i++) { + psVector *values = psVectorAlloc(0, PS_TYPE_F32); + pixels->data[i] = values; + } + + // Pull the pixels out into the vectors + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) { + // the overscan and image might not be aligned. pixels->data represents + // the image row pixels. + int diff = overscan->col0 - image->col0; // Offset between the two regions + for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) { + int iFixed = i - diff; + // i is column on image + // iFixed is column on overscan + psVector *values = pixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + overscan->numRows); + for (int j = 0; j < overscan->numRows; j++) { + values->data.F32[index++] = overscan->data.F32[j][iFixed]; + } + values->n += overscan->numRows; + pixels->data[i] = values; // Update the pointer in case it's moved + } + } + psFree(iter); + + // Reduce the overscans + psVector *reduced = pmOverscanVector(&chi2, overscanOpts, pixels, stats); + psFree(pixels); + if (! reduced) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); + return false; + } + + // Subtract column by column + for (int i = 0; i < image->numCols; i++) { + for (int j = 0; j < image->numRows; j++) { + image->data.F32[j][i] -= reduced->data.F32[i]; + } + } + psFree(reduced); + } + + pmOverscanUpdateHeader (hdu, overscanOpts, chi2); + psFree(stats); + + return true; + +} // End of overscan subtraction + Index: /trunk/psModules/src/detrend/pmOverscan.h =================================================================== --- /trunk/psModules/src/detrend/pmOverscan.h (revision 14505) +++ /trunk/psModules/src/detrend/pmOverscan.h (revision 14505) @@ -0,0 +1,72 @@ +/* @file pmOverscan.h + * @brief Functions to subtract the overscan, used by pmBiasSubtract + * + * @author George Gusciora, MHPCC + * @author Paul Price, IfA + * @author Eugene Magnier, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-15 20:21:18 $ + * Copyright 2004--2006 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_OVERSCAN_H +#define PM_OVERSCAN_H + +/// @addtogroup detrend Detrend Creation and Application +/// @{ + +/// Type of fit to perform +typedef enum { + PM_FIT_NONE, ///< No fit + PM_FIT_POLY_ORD, ///< Fit ordinary polynomial + PM_FIT_POLY_CHEBY, ///< Fit Chebyshev polynomial + PM_FIT_SPLINE ///< Fit cubic splines +} pmFit; + +/// Options for overscan subtraction +/// +/// The overscan subtraction may be performed by reducing all overscan regions to a single value (e.g., if +/// there is no structure); or the overscan may be fit perpendicular to the read direction (usually the +/// columns) with a particular functional form; or a single value may be subtracted for each read/scan without +/// fitting (if the structure defies characterisation). In any case, statistics are required to reduce +/// multiple values to a single value (either for the scan, or for the entire overscan regions). +typedef struct +{ + // Inputs + bool single; ///< Reduce all overscan regions to a single value? + bool constant; ///< use a supplied constant value (do not measure region) + float value; ///< supplied value if needed (per above) + pmFit fitType; ///< Type of fit to overscan + unsigned int order; ///< Order of polynomial, or number of spline pieces + psStats *stat; ///< Statistic to use when reducing the minor direction + int boxcar; ///< Boxcar smoothing radius + float gauss; ///< Gaussian smoothing sigma + // Outputs + psPolynomial1D *poly; ///< Result of polynomial fit + psSpline1D *spline; ///< Result of spline fit +} +pmOverscanOptions; + +/// Allocator for overscan options +pmOverscanOptions *pmOverscanOptionsAlloc(bool single, ///< Reduce all overscan regions to a single value? + pmFit fitType, ///< Type of fit to overscan + unsigned int order, ///< Order of polynomial, or number of splines + psStats *stat, ///< Statistic to use + int boxcar, ///< Boxcar smoothing radius + float gauss ///< Gaussian smoothing sigma + ); + +psVector *pmOverscanVector(float *chi2, // chi^2 from fit + pmOverscanOptions *overscanOpts, // Overscan options + const psArray *pixels, // Array of vectors containing the pixel values + psStats *myStats // Statistic to use in reducing the overscan + ); + +bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2); + +bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts); + +/// @} +#endif + Index: /trunk/psModules/src/objects/pmPSFtry.h =================================================================== --- /trunk/psModules/src/objects/pmPSFtry.h (revision 14504) +++ /trunk/psModules/src/objects/pmPSFtry.h (revision 14505) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-06-20 02:22:26 $ + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-15 20:21:18 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -107,6 +107,5 @@ */ bool pmPSFtryMetric_Alt( - pmPSFtry *try - , ///< Add comment. + pmPSFtry *try, ///< Add comment. float RADIUS ///< Add comment. ); Index: /trunk/psModules/src/objects/pmSource.h =================================================================== --- /trunk/psModules/src/objects/pmSource.h (revision 14504) +++ /trunk/psModules/src/objects/pmSource.h (revision 14505) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-07-20 20:12:59 $ + * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-15 20:21:18 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -77,8 +77,8 @@ float extMag; ///< calculated from flux in modelEXT float errMag; ///< error in psfMag OR extMag (depending on type) - float apMag; ///< apMag corresponding to psfMag or extMag (depending on type) - float pixWeight; // model-weighted coverage of valid pixels - psRegion region; // area on image covered by selected pixels - float sky, skyErr; //?< The sky and its error at the center of the object + float apMag; ///< apMag corresponding to psfMag or extMag (depending on type) + float pixWeight; ///< model-weighted coverage of valid pixels + psRegion region; ///< area on image covered by selected pixels + float sky, skyErr; ///< The sky and its error at the center of the object } pmSource; Index: /trunk/psModules/src/psmodules.h =================================================================== --- /trunk/psModules/src/psmodules.h (revision 14504) +++ /trunk/psModules/src/psmodules.h (revision 14505) @@ -61,4 +61,5 @@ #include #include +#include #include #include