Index: trunk/psModules/src/detrend/Makefile.am =================================================================== --- trunk/psModules/src/detrend/Makefile.am (revision 5675) +++ trunk/psModules/src/detrend/Makefile.am (revision 6872) @@ -3,15 +3,18 @@ libpsmoduledetrend_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS) libpsmoduledetrend_la_LDFLAGS = -release $(PACKAGE_VERSION) -libpsmoduledetrend_la_SOURCES = pmFlatField.c \ - pmMaskBadPixels.c \ - pmNonLinear.c +libpsmoduledetrend_la_SOURCES = \ + pmFlatField.c \ + pmFringeStats.c \ + pmMaskBadPixels.c \ + pmNonLinear.c psmoduleincludedir = $(includedir) psmoduleinclude_HEADERS = \ - pmFlatField.h \ - pmFlatFieldErrors.h \ - pmMaskBadPixelsErrors.h \ - pmMaskBadPixels.h \ - pmNonLinear.h + pmFlatField.h \ + pmFlatFieldErrors.h \ + pmFringeStats.h \ + pmMaskBadPixelsErrors.h \ + pmMaskBadPixels.h \ + pmNonLinear.h EXTRA_DIST = pmFlatFieldErrors.dat pmMaskBadPixelsErrors.dat Index: trunk/psModules/src/detrend/pmFlatField.c =================================================================== --- trunk/psModules/src/detrend/pmFlatField.c (revision 5675) +++ trunk/psModules/src/detrend/pmFlatField.c (revision 6872) @@ -1,2 +1,8 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + /** @file pmFlatField.c * @@ -18,6 +24,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-18 19:43:14 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -30,5 +36,5 @@ #include #include -#include +#include #include "pslib.h" @@ -36,32 +42,8 @@ #include "pmMaskBadPixels.h" #include "pmFlatFieldErrors.h" -#include "pmSubtractBias.h" -// XXX: This should be removed when the autoconf stuff handles psConstants.h correctly. -#define PS_WARN_PTR_NON_NULL(NAME) \ -if ((NAME) == NULL) { \ - psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \ -} \ -bool pmFlatField( - pmReadout *in, - const pmReadout *flat) +bool pmFlatField(pmReadout *in, const pmReadout *flat) { - // XXX: Use the proper image and readout asserts. - PS_ASSERT_PTR_NON_NULL(in, true); - PS_ASSERT_PTR_NON_NULL(in->image, false); - PS_ASSERT_PTR_NON_NULL(in->mask, false); - PS_ASSERT_PTR_NON_NULL(flat, false); - PS_ASSERT_PTR_NON_NULL(flat->image, false); - if (in == NULL) - printf("XXX: NULL\n"); - - // XXX: Not sure if this is correct. Must consult with IfA. - PS_ASSERT_PTR_NON_NULL(in->mask, false); - - PS_WARN_PTR_NON_NULL(in->parent); - if (in->parent != NULL) { - PS_WARN_PTR_NON_NULL(in->parent->concepts); - } int i = 0; int j = 0; @@ -71,28 +53,43 @@ psElemType flatType; psElemType maskType; - psImage *inMask = NULL; - psImage *flatImage = NULL; - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(in); - flatImage = flat->image; - inMask = in->mask; + // Check for nulls + if (in == NULL) { + return true; // Readout may not have data in it + } else if(flat==NULL) { + psError( PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_pmFlatField_NULL_FLAT_READOUT); + return false; + } + + psImage *inImage = in->image; // Input image + psImage *inMask = in->mask; // Mask for input image + psImage *flatImage = flat->image; // Flat-field image + + // Offsets on the chip + int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0"); + int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0"); + int x0flat = psMetadataLookupS32(NULL, flat->parent->concepts, "CELL.X0"); + int y0flat = psMetadataLookupS32(NULL, flat->parent->concepts, "CELL.Y0"); + + if (inImage == NULL) { + psError( PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_pmFlatField_NULL_INPUT_IMAGE); + return false; + } else if(flatImage == NULL) { + psError( PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_pmFlatField_NULL_FLAT_IMAGE); + return false; + } // Check input image and its mask are not larger than flat image - if (trimmedImg == NULL) - printf("XXX: 00\n"); - if (flatImage == NULL) - printf("XXX 01\n"); - - if (trimmedImg->numRows>flatImage->numRows || trimmedImg->numCols>flatImage->numCols) { + if (inImage->numRows>flatImage->numRows || inImage->numCols>flatImage->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE, - trimmedImg->numRows, trimmedImg->numCols, flatImage->numRows, flatImage->numCols); + inImage->numRows, inImage->numCols, flatImage->numRows, flatImage->numCols); return false; } - if (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols) { + if (inMask && (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols)) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmFlatField_SIZE_MASK_IMAGE, @@ -102,6 +99,6 @@ // Determine total offset based on image offset with chip offset - totOffCol = trimmedImg->col0 + in->col0; - totOffRow = trimmedImg->row0 + in->row0; + totOffCol = inImage->col0 + y0in - flatImage->col0 - y0flat; + totOffRow = inImage->row0 + x0in - flatImage->row0 - x0flat; // Check that offsets are within image limits @@ -111,10 +108,10 @@ totOffRow, totOffCol, flatImage->numRows, flatImage->numCols); return false; - } else if(totOffRow>=trimmedImg->numRows || totOffCol>=trimmedImg->numCols) { + } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE, - totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols); + totOffRow, totOffCol, inImage->numRows, inImage->numCols); return false; - } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) { + } else if(inMask && (totOffRow>=inMask->numRows || totOffCol>=inMask->numCols)) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmFlatField_OFFSET_MASK_IMAGE, @@ -124,5 +121,5 @@ // Check for incorrect types - inType = trimmedImg->type.type; + inType = inImage->type.type; flatType = flatImage->type.type; maskType = inMask->type.type; @@ -137,5 +134,5 @@ flatType); return false; - } else if(maskType != PS_TYPE_MASK) { + } else if(inMask && inMask->type.type != PS_TYPE_MASK) { psError( PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_pmFlatField_TYPE_MASK_IMAGE, @@ -153,18 +150,20 @@ case PS_TYPE_##TYPE: \ /* Per Eugene's request, use two sets of loops: first to fill mask, second to avoid div with bad pix */ \ - for(j = totOffRow; j < trimmedImg->numRows; j++) { \ - for(i = totOffCol; i < trimmedImg->numCols; i++) { \ + for(j = totOffRow; j < inImage->numRows; j++) { \ + for(i = totOffCol; i < inImage->numCols; i++) { \ if(flatImage->data.TYPE[j][i] <= 0.0) { \ /* Negative or zero flat pixels shall be masked in input image as PM_MASK_FLAT */ \ - inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT; \ + if (inMask) { \ + inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT; \ + } \ flatImage->data.TYPE[j][i] = 0.0; \ } \ } \ } \ - for(j = totOffRow; j < trimmedImg->numRows; j++) { \ - for(i = totOffCol; i < trimmedImg->numCols; i++) { \ - if(!inMask->data.PS_TYPE_MASK_DATA[j][i]) { \ + for(j = totOffRow; j < inImage->numRows; j++) { \ + for(i = totOffCol; i < inImage->numCols; i++) { \ + if(inMask && !inMask->data.PS_TYPE_MASK_DATA[j][i]) { \ /* Module shall divide the input image by the flat-fielded image */ \ - trimmedImg->data.TYPE[j][i] /= flatImage->data.TYPE[j][i]; \ + inImage->data.TYPE[j][i] /= flatImage->data.TYPE[j][i]; \ } \ } \ Index: trunk/psModules/src/detrend/pmFlatField.h =================================================================== --- trunk/psModules/src/detrend/pmFlatField.h (revision 5675) +++ trunk/psModules/src/detrend/pmFlatField.h (revision 6872) @@ -1,2 +1,8 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + /** @file pmFlatField.h * @@ -18,6 +24,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-20 23:06:24 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -25,5 +31,5 @@ #include "pslib.h" -#include "pmAstrometry.h" +#include "pmFPA.h" Index: trunk/psModules/src/detrend/pmFringeStats.c =================================================================== --- trunk/psModules/src/detrend/pmFringeStats.c (revision 6872) +++ trunk/psModules/src/detrend/pmFringeStats.c (revision 6872) @@ -0,0 +1,240 @@ +/** @file pmFringeStats.c + * + * @author Eugene Magnier, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 IfA + */ + +#include "pslib.h" +#include "pmFPA.h" +#include "pmFringeStats.h" + +static void fringeRegionsFree(pmFringeRegions *fringe) +{ + psFree(fringe->x); + psFree(fringe->y); + psFree(fringe->mask); + return; +} + +pmFringeRegions *pmFringeRegionsAlloc(int nPts, int dX, int dY, int nX, int nY) +{ + pmFringeRegions *fringe = psAlloc(sizeof(pmFringeRegions)); + (void)psMemSetDeallocator(fringe, (psFreeFunc)fringeRegionsFree); + + fringe->x = NULL; + fringe->y = NULL; + fringe->mask = NULL; + + fringe->nRequested = nPts; + fringe->nAccepted = 0; + + fringe->dX = dX; + fringe->dY = dY; + fringe->nX = nX; + fringe->nY = nY; + + return fringe; +} + + +static void fringeStatsFree(pmFringeStats *stats) +{ + psFree(stats->regions); + psFree(stats->f); + psFree(stats->df); +} + +pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions) +{ + pmFringeStats *stats = psAlloc(sizeof(pmFringeStats)); + (void)psMemSetDeallocator(stats, (psFreeFunc)fringeStatsFree); + + stats->regions = psMemIncrRefCounter(regions); + stats->f = psVectorAlloc(regions->x->n, PS_TYPE_F32); + stats->df = psVectorAlloc(regions->x->n, PS_TYPE_F32); + + return stats; +} + + +static void fringeScaleFree(pmFringeScale *scale) +{ + psFree(scale->coeff); + psFree(scale->coeffErr); + return; +} + +pmFringeScale *pmFringeScaleAlloc(int nFringeFrames) +{ + pmFringeScale *scale = psAlloc(sizeof(pmFringeScale)); + (void)psMemSetDeallocator(scale, (psFreeFunc)fringeScaleFree); + + scale->nFringeFrames = nFringeFrames; + scale->coeff = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32); + scale->coeffErr = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32); + + return scale; +} + +bool pmFringeStatsCreatePoints(pmFringeRegions *fringe, psImage *image) +{ + + double frnd; + // create fringe->nRequested + + psRandom *rnd = psRandomAlloc(PS_RANDOM_TAUS, 0); + + fringe->x = psVectorRecycle(fringe->x, fringe->nRequested, PS_TYPE_F32); + fringe->y = psVectorRecycle(fringe->y, fringe->nRequested, PS_TYPE_F32); + fringe->mask = psVectorRecycle(fringe->mask, fringe->nRequested, PS_TYPE_U8); + + int nX = image->numCols; + int nY = image->numRows; + + psF32 *xPt = fringe->x->data.F32; + psF32 *yPt = fringe->y->data.F32; + + int dX = fringe->dX; + int dY = fringe->dY; + + // generate random points located within image bounds + for (int i = 0; i < fringe->nRequested; i++) { + frnd = psRandomUniform(rnd); + xPt[i] = (nX - 2*dX)* frnd + dX; + frnd = psRandomUniform(rnd); + yPt[i] = (nY - 2*dY)* frnd + dY; + } + return true; +} + +pmFringeStats *pmFringeStatsMeasure(pmFringeRegions *fringe, pmReadout *readout, psMaskType maskVal) +{ + if (!fringe->x || !fringe->y) { + // create the fringe vectors for this image + pmFringeStatsCreatePoints(fringe, readout->image); + } + + PS_ASSERT_PTR_NON_NULL(fringe->x, false); + PS_ASSERT_PTR_NON_NULL(fringe->y, false); + + pmFringeStats *measurements = pmFringeStatsAlloc(fringe); + + psF32 *xPt = fringe->x->data.F32; + psF32 *yPt = fringe->y->data.F32; + psF32 *fPt = measurements->f->data.F32; + psF32 *dfPt = measurements->df->data.F32; + + int dX = fringe->dX; + int dY = fringe->dY; + + psImage *image = readout->image; + psImage *mask = readout->mask; + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + + for (int i = 0; i < fringe->x->n; i++) { + psRegion region = psRegionSet(xPt[i] - dX, yPt[i] - dY, xPt[i] + dX + 1, yPt[i] + dY + 1); + region = psRegionForImage(image, region); + psImage *subImage = psImageSubset(image, region); + psImage *subMask = psImageSubset(mask, region); + psImageStats(stats, subImage, subMask, maskVal); + + fPt[i] = stats->sampleMedian; + dfPt[i] = stats->sampleStdev; + } + + return measurements; +} + +// XXX include the fringe error (fringe->df) in the fit? +pmFringeScale *pmFringeScaleMeasure(pmFringeStats *science, psArray *fringes) +{ + double sum; + psVector *v1; + psVector *v2; + + pmFringeScale *scale = pmFringeScaleAlloc(fringes->n); + + int nCof = fringes->n + 1; + int nPts = science->f->n; + + psImage *A = psImageAlloc(nCof, nCof, PS_TYPE_F64); + psVector *B = psVectorAlloc(nCof, PS_TYPE_F64); + + for (int i = 0; i < nCof; i++) { + psF32 *p1 = NULL; + if (i != 0) { + v1 = fringes->data[i - 1]; + p1 = v1->data.F32; + } + for (int j = 0; j < nCof; j++) { + psF32 *p2 = NULL; + if (j != 0) { + v2 = fringes->data[j - 1]; + p2 = v2->data.F32; + } + + sum = 0; + for (int k = 0; k < nPts; k++) { + psF32 f1 = (p1 == NULL) ? 1.0 : p1[k]; + psF32 f2 = (p2 == NULL) ? 1.0 : p2[k]; + sum += f1*f2; + } + A->data.F32[i][j] = sum; + } + + sum = 0; + psF32 *p2 = science->f->data.F32; + for (int j = 0; j < nPts; j++) { + psF32 f1 = (p1 == NULL) ? 1.0 : p1[j]; + psF32 f2 = p2[j]; + sum += f1*f2; + } + B->data.F32[i] = sum; + } + + if (! psGaussJordan(A, B)) { + psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations. Returning NULL.\n"); + } + + + for (int i = 0; i < nCof; i++) { + scale->coeff->data.F32[i] = B->data.F32[i]; + scale->coeffErr->data.F32[i] = sqrt(A->data.F32[i][i]); + } + + return scale; +} + +// XXX note that this modifies the input fringe images +psImage *pmFringeCorrect(pmReadout *readout, pmFringeRegions *fringes, psArray *fringeImages, psArray *fringeStats, psMaskType maskVal) +{ + // measure the fringe stats for the science frame and solve for the scales + pmFringeStats *scienceStats = pmFringeStatsMeasure(fringes, readout, maskVal); + pmFringeScale *scale = pmFringeScaleMeasure(scienceStats, fringeStats); + psFree(scienceStats); + + // build the fringe correction image + // XXX we could save data space by making the first image the output image + psImage *sumFringe = NULL; + for (int i = 0; i < fringeImages->n; i++) { + + // rescale the fringe image + psBinaryOp(fringeImages->data[i], fringeImages->data[i], "*", + psScalarAlloc(scale->coeff->data.F32[i+1], PS_TYPE_F32)); + + // sum together + sumFringe = (psImage*)psBinaryOp(sumFringe, sumFringe, "+", fringeImages->data[i]); + } + psFree(scale); + + // subtract the resulting fringe frame + readout->image = (psImage*)psBinaryOp(readout->image, readout->image, "-", sumFringe); + + return sumFringe; +} + Index: trunk/psModules/src/detrend/pmFringeStats.h =================================================================== --- trunk/psModules/src/detrend/pmFringeStats.h (revision 6872) +++ trunk/psModules/src/detrend/pmFringeStats.h (revision 6872) @@ -0,0 +1,92 @@ +/** @file pmFringeStats.h + * + * @brief Measure Fringe statistics and apply correction + * + * @author Eugene Magnier, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 IfA, University of Hawaii + */ + +# ifndef PM_FRINGE_STATS +# define PM_FRINGE_STATS + +/** Structure to hold the fringe measurement regions. + */ +typedef struct +{ + psU32 nRequested; // number of fringe points selected + psU32 nAccepted; // number of fringe points not masked + psU32 dX; // median box half-width + psU32 dY; // median box half-height + psU32 nX; // large-scale smoothing in x (col) + psU32 nY; // large-scale smoothing in y (row) + psVector *x; // fringe point coordinates (col) + psVector *y; // fringe point coordinates (row) + psVector *mask; // fringe point on/off mask +} +pmFringeRegions; + +pmFringeRegions *pmFringeRegionsAlloc ( + int nPts, // number of points to create + int dX, // half-width of fringe boxes + int dY, // half-height of fringe boxes + int nX, // smoothing scale in x + int nY // smoothing scale in y +); + +/** Structure to hold the fringe measurements for a particular image + */ +typedef struct +{ + pmFringeRegions *regions; // Fringe regions + psVector *f; // fringe point median + psVector *df; // fringe point stdev +} +pmFringeStats; + +pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions); + + +/** the pmFringeScale structure defines the relationship between two fringe measurements + */ +typedef struct +{ + int nFringeFrames; + psVector *coeff; + psVector *coeffErr; +} +pmFringeScale; + +/** Measure the fringe stats for an image + * + * Given an input image and a fringe stats structure, measure the statistics for each of the fringe points on + * the image. If the fringe vectors are NULL, new vectors are created. + * + * XXX should the structure carry the image dimensions for validation? + * + * @return bool: True or false for success or failure + */ +pmFringeStats *pmFringeStatsMeasure( + pmFringeRegions *fringe, ///< fringe regions + pmReadout *readout, ///< measure fringes on this readout + psMaskType maskVal ///< Mask value for statistics +); + +/** Fringe correct the science image + * + * Given a science image, a set of fringe stats and a matched set of fringe images, correct the science image + * for the fringe images. + * + * @return psImage: summed fringe image + */ +psImage *pmFringeCorrect(pmReadout *in, // place info about results here + pmFringeRegions *fringes, // The fringe regions used + psArray *fringeImages, // fringe images to use in correction + psArray *fringeStats, // fringe stats to use in correction + psMaskType maskVal // Value to mask + ); + +# endif Index: trunk/psModules/src/detrend/pmMaskBadPixels.c =================================================================== --- trunk/psModules/src/detrend/pmMaskBadPixels.c (revision 5675) +++ trunk/psModules/src/detrend/pmMaskBadPixels.c (revision 6872) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmMaskBadPixels.c * @@ -19,6 +24,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-18 19:43:14 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -31,157 +36,77 @@ #include #include -#include +#include #include "pmMaskBadPixels.h" #include "pmMaskBadPixelsErrors.h" -#include "pmSubtractBias.h" -//XXX: REmove, autoconf is broken. -#define PS_WARN_PTR_NON_NULL(NAME) \ -if ((NAME) == NULL) { \ - psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \ -} \ +bool pmMaskBadPixels(pmReadout *in, const psImage *mask, unsigned int maskVal, float sat, + unsigned int growVal, int grow) +{ + int i = 0; + int j = 0; + int jj = 0; + int ii = 0; + int rRound = 0; + int rowMin = 0; + int rowMax = 0; + int colMin = 0; + int colMax = 0; + int totOffCol = 0; + int totOffRow = 0; + float r = 0.0f; + psElemType inType; + psElemType maskType; + psImage *inImage = NULL; + psImage *inMask = NULL; -/****************************************************************************** -GrowPixel(inMask, row, col, radius, maskVal): This private routine takes an -input image mask and a pixel location, then sets (logical or) all pixels with -parameter radius if that pixel to maskVal parameter. -*****************************************************************************/ -psBool GrowPixel( - psImage *inMask, - psS32 col, - psS32 row, - psS32 radius, - psU32 maskVal) -{ - psS32 rowMin = PS_MAX(row-radius, 0); - psS32 rowMax = PS_MIN(row+radius+1, inMask->numRows); - psS32 colMin = PS_MAX(col-radius, 0); - psS32 colMax = PS_MIN(col+radius+1, inMask->numCols); - psF32 squareRadius = PS_SQR(radius); + // Check for nulls + if (in == NULL) { + return true; // Readout may not have data in it + } else if(mask==NULL) { + psError( PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_pmMaskBadPixels_NULL_MASK_IMAGE); + return false; + } - - for(psS32 i=rowMin; idata.U8[i][j] |= maskVal; - } - } - } - return(true); -} - -/****************************************************************************** -GetRadius(inImg, col, row, sat, growVal, grow): This private routine takes an -input image and pixel location and determines what radius that pixel must grow -by. - -//XXX: Inline this or macro it for speed. - *****************************************************************************/ -static psS32 GetRadius( - psImage *inImg, - psS32 col, - psS32 row, - psF32 sat, - psU32 growVal, - psS32 grow) -{ - psS32 growRadius = 0; - if (inImg->type.type == PS_TYPE_F32) { - if(inImg->data.F32[row][col] == growVal) { - growRadius = grow; - } - if (inImg->data.F32[row][col] > sat) { - growRadius = PS_MAX(growRadius+1, 2); - } - } else if (inImg->type.type == PS_TYPE_S32) { - if(inImg->data.S32[row][col] == growVal) { - growRadius = grow; - } - if (inImg->data.S32[row][col] > sat) { - growRadius = PS_MAX(growRadius+1, 2); - } - } else if (inImg->type.type == PS_TYPE_U16) { - if(inImg->data.U16[row][col] == growVal) { - growRadius = grow; - } - if (inImg->data.U16[row][col] > sat) { - growRadius = PS_MAX(growRadius+1, 2); - } - } else { - psError( PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_pmMaskBadPixels_TYPE_MASK_IMAGE, - inImg->type.type); - } - return(growRadius); -} - - -bool pmMaskBadPixels( - pmReadout *in, - const pmReadout *mask, - unsigned int maskVal, - float sat, - unsigned int growVal, - int grow) -{ - // XXX: Review this code then put all asserts at the top. - PS_ASSERT_PTR_NON_NULL(in, true); - PS_ASSERT_PTR_NON_NULL(in->image, false); - PS_WARN_PTR_NON_NULL(in->parent); - if (in->parent != NULL) { - PS_WARN_PTR_NON_NULL(in->parent->concepts); - } - PS_ASSERT_PTR_NON_NULL(mask, false); - int i = 0; - int j = 0; - int totOffCol = 0; - int totOffRow = 0; - psElemType inType; - psElemType maskType; - psImage *inMask = NULL; - - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(in); - if(in->mask == NULL) { - in->mask = psImageAlloc(trimmedImg->numCols, trimmedImg->numRows, PS_TYPE_MASK); - PS_IMAGE_SET_U8(in->mask, 0); + inImage = in->image; + if (inImage == NULL) { + psError( PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_pmMaskBadPixels_NULL_INPUT_IMAGE); + return false; + } else if(in->mask == NULL) { + in->mask = psImageAlloc(inImage->numCols, inImage->numRows, PS_TYPE_MASK); + memset(in->mask->data.V[0], 0, inImage->numCols*inImage->numRows*PSELEMTYPE_SIZEOF(PS_TYPE_MASK)); } inMask = in->mask; // Check input image and its mask are not larger than mask - if(trimmedImg->numRows > mask->image->numRows || trimmedImg->numCols > mask->image->numCols) { + if(inImage->numRows > mask->numRows || inImage->numCols > mask->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmMaskBadPixels_SIZE_INPUT_IMAGE, - trimmedImg->numRows, trimmedImg->numCols, mask->image->numRows, mask->image->numCols); + inImage->numRows, inImage->numCols, mask->numRows, mask->numCols); return false; - } else if(inMask->numRows > mask->image->numRows || inMask->numCols > mask->image->numCols) { + } else if(inMask->numRows>mask->numRows || inMask->numCols>mask->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmMaskBadPixels_SIZE_MASK_IMAGE, - inMask->numRows, inMask->numCols, mask->image->numRows, mask->image->numCols); + inMask->numRows, inMask->numCols, mask->numRows, mask->numCols); return false; } // Determine total offset based on image offset with chip offset - totOffCol = trimmedImg->col0 + in->col0; - totOffRow = trimmedImg->row0 + in->row0; + totOffCol = inImage->col0; // + in->col0; + totOffRow = inImage->row0; // + in->row0; // Check that offsets are within image limits - if(totOffRow>=mask->image->numRows || totOffCol>=mask->image->numCols) { + if(totOffRow>=mask->numRows || totOffCol>=mask->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmMaskBadPixels_OFFSET_MASK_IMAGE, - totOffRow, totOffCol, mask->image->numRows, mask->image->numCols); + totOffRow, totOffCol, mask->numRows, mask->numCols); return false; - } else if(totOffRow>=trimmedImg->numRows || totOffCol>=trimmedImg->numCols) { + } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) { psError( PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE, - totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols); + totOffRow, totOffCol, inImage->numRows, inImage->numCols); return false; } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) { @@ -193,6 +118,6 @@ // Check for incorrect types - inType = trimmedImg->type.type; - maskType = mask->image->type.type; + inType = inImage->type.type; + maskType = mask->type.type; if(PS_IS_PSELEMTYPE_COMPLEX(inType)) { psError( PS_ERR_BAD_PARAMETER_TYPE, true, @@ -207,34 +132,57 @@ } - // - // This is the main loop which examines each pixel and masks pixels if - // A: The mask matches the maskValue - // B: The pixel is larger than the saturation value - // C: The pixel equals the grow value (in which case a circle is masked) - // - for(i=totOffRow; inumRows; i++) { - for(j=totOffCol; jnumCols; j++) { - // - // A: Pixels which satisfy maskVal shall be masked - // - if (mask->image->data.U8[i][j] == maskVal) { - in->mask->data.U8[i][j] |= maskVal; - } + // Macro for all PS types + #define PM_BAD_PIXELS(TYPE) \ +case PS_TYPE_##TYPE: \ + for(j=totOffRow; jnumRows; j++) { \ + for(i=totOffCol; inumCols; i++) { \ + \ + /* Pixels with flux greater than sat shall be masked */ \ + if(inImage->data.TYPE[j][i] > sat) { \ + inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_SAT; \ + } \ + \ + /* Pixels which satisfy maskVal shall be masked */ \ + inMask->data.PS_TYPE_MASK_DATA[j][i] |= (mask->data.PS_TYPE_MASK_DATA[j][i]&maskVal); \ + \ + /* Pixels which satisfy growVal and within the grow radius shall be masked */ \ + if(mask->data.PS_TYPE_MASK_DATA[j][i] & growVal) { \ + rowMin = MAX(j-grow, 0); \ + rowMax = MIN(j+grow+1, inImage->numRows); \ + colMin = MAX(i-grow, 0); \ + colMax = MIN(i+grow+1, inImage->numCols); \ + for(jj=rowMin; jjdata.PS_TYPE_MASK_DATA[jj][ii] |= \ + (mask->data.PS_TYPE_MASK_DATA[j][i]&growVal); \ + } \ + } \ + } \ + } \ + } \ + } \ + break; - // - // We first determine how much we need to grow by and store this in - // growRadius. - // - psS32 growRadius = GetRadius(trimmedImg, j, i, sat, growVal, grow); - - // - // Grow the pixel mask if necessary. - // - if (growRadius != 0) { - GrowPixel(in->mask, j, i, growRadius, maskVal); - } - } + // Switch to call bad pixel masking macro defined above + switch(inType) { + PM_BAD_PIXELS(U8); + PM_BAD_PIXELS(U16); + PM_BAD_PIXELS(U32); + PM_BAD_PIXELS(U64); + PM_BAD_PIXELS(S8); + PM_BAD_PIXELS(S16); + PM_BAD_PIXELS(S32); + PM_BAD_PIXELS(S64); + PM_BAD_PIXELS(F32); + PM_BAD_PIXELS(F64); + default: + psError( PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_pmMaskBadPixels_TYPE_UNSUPPORTED, + inType); } - return true; + return false; } Index: trunk/psModules/src/detrend/pmMaskBadPixels.h =================================================================== --- trunk/psModules/src/detrend/pmMaskBadPixels.h (revision 5675) +++ trunk/psModules/src/detrend/pmMaskBadPixels.h (revision 6872) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmMaskBadPixels.h * @@ -19,6 +24,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-15 20:09:03 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -26,5 +31,5 @@ #include "pslib.h" -#include "pmAstrometry.h" +#include "pmFPA.h" /** Mask values */ @@ -33,8 +38,8 @@ PM_MASK_BADCOL = 0x0002, ///< The pixel is a bad column. PM_MASK_SAT = 0x0004, ///< The pixel is saturated. - PM_MASK_FLAT = 0x0008 ///< The pixel is non-positive in the flat-field. + PM_MASK_BAD = 0x0008, ///< The pixel is low + PM_MASK_FLAT = 0x0010 ///< The pixel is non-positive in the flat-field. } pmMaskValue; -// XXX: Use PS_MIN, PS_MAX /** Macro to find maximum of two numbers */ #define MAX(A,B)((A)>=(B)?(A):(B)) @@ -54,5 +59,5 @@ bool pmMaskBadPixels( pmReadout *in, ///< Readout containing input image data. - const pmReadout *mask, ///< Mask data to be added to readout mask data. + const psImage *mask, ///< Mask data to be added to readout mask data. unsigned int maskVal, ///< Mask value to determine what to add to input mask. float sat, ///< Saturation limit to mask bad pixels. Index: trunk/psModules/src/detrend/pmMaskBadPixelsErrors.h =================================================================== --- trunk/psModules/src/detrend/pmMaskBadPixelsErrors.h (revision 5675) +++ trunk/psModules/src/detrend/pmMaskBadPixelsErrors.h (revision 6872) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmMaskBadPixelsErrors.h * @@ -7,6 +12,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-28 20:43:52 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -23,5 +28,5 @@ * $2 The error text (rest of the line in pmMaskBadPixelsErrors.h) * $n The order of the source line in pmMaskBadPixelsErrors.h (comments excluded) - * + * * DO NOT EDIT THE LINES BETWEEN //~Start and //~End! ANY CHANGES WILL BE OVERWRITTEN. */ Index: trunk/psModules/src/detrend/pmNonLinear.c =================================================================== --- trunk/psModules/src/detrend/pmNonLinear.c (revision 5675) +++ trunk/psModules/src/detrend/pmNonLinear.c (revision 6872) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmNonLinear.c * @@ -5,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-12-05 20:49:40 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -22,22 +27,14 @@ #include "pmNonLinear.h" -#include "pmSubtractBias.h" -// XXX: Remove, autoconf must be -#define PS_WARN_PTR_NON_NULL(NAME) \ -if ((NAME) == NULL) { \ - psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \ -} \ /****************************************************************************** pmNonLinearityLookup(): This routine will take an pmReadout image as input and a 1-D polynomial. For each pixel in the input image, the polynomial will -be evaluated at that pixels value, and the image pixel will then be set -to +be evaluated at that pixels value, and the image pixel will then be set to that value. -*****************************************************************************/ + *****************************************************************************/ -pmReadout *pmNonLinearityPolynomial( - pmReadout *inputReadout, - const psPolynomial1D *input1DPoly) +pmReadout *pmNonLinearityPolynomial(pmReadout *inputReadout, + const psPolynomial1D *input1DPoly) { PS_ASSERT_PTR_NON_NULL(inputReadout, NULL); @@ -45,18 +42,11 @@ PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL); PS_ASSERT_PTR_NON_NULL(input1DPoly, NULL); - PS_WARN_PTR_NON_NULL(inputReadout->parent); - if (inputReadout->parent != NULL) { - PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts); - } - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout); + psS32 i; + psS32 j; - for (psS32 i=0;inumRows;i++) { - for (psS32 j=0;jnumCols;j++) { - trimmedImg->data.F32[i][j] = psPolynomial1DEval(input1DPoly, - trimmedImg->data.F32[i][j]); + for (i=0;iimage->numRows;i++) { + for (j=0;jimage->numCols;j++) { + inputReadout->image->data.F32[i][j] = psPolynomial1DEval(input1DPoly, inputReadout->image->data.F32[i][j]); } } @@ -71,48 +61,68 @@ inFluxe, and the corresponding value in outFlux. The image pixel will then be set to the value from outFlux. - -XXX: Must assert that filename exists. This should probably happen in -the lookup files. *****************************************************************************/ -pmReadout *pmNonLinearityLookup( - pmReadout *inputReadout, - const char *filename -) +pmReadout *pmNonLinearityLookup(pmReadout *inputReadout, + const psVector *inFlux, + const psVector *outFlux) { PS_ASSERT_PTR_NON_NULL(inputReadout,NULL); PS_ASSERT_PTR_NON_NULL(inputReadout->image,NULL); PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL); - PS_WARN_PTR_NON_NULL(inputReadout->parent); - if (inputReadout->parent != NULL) { - PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts); + PS_ASSERT_PTR_NON_NULL(inFlux,NULL); + if (inFlux->n < 2) { + psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): input vector less than 2 elements. Returning inputReadout image."); + return(inputReadout); } - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout); + PS_ASSERT_PTR_NON_NULL(outFlux,NULL); + psS32 tableSize = inFlux->n; + if (inFlux->n != outFlux->n) { + tableSize = PS_MIN(inFlux->n, outFlux->n); + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmNonLinear.c: pmNonLinearityLookup(): input vectors have different sizes (%d, %d)\n", inFlux->n, outFlux->n); + } + PS_ASSERT_VECTOR_TYPE(inFlux, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(outFlux, PS_TYPE_F32, NULL); - psLookupTable *tmpLT = psLookupTableAlloc(filename, "%f %f", 0); - psS32 numLines = psLookupTableRead(tmpLT); + psS32 i; + psS32 j; + psS32 binNum; + psScalar x; psS32 numPixels = 0; - if (numLines < 2) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: Lookup Table is too small. Returning original pmReadout.\n"); - } else { - for (psS32 i=0;inumRows;i++) { - for (psS32 j=0;jnumCols;j++) { - psF64 tmpD = psLookupTableInterpolate(tmpLT, trimmedImg->data.F32[i][j], 1); - if (!isnan(tmpD)) { - trimmedImg->data.F32[i][j] = tmpD; - } else { - numPixels++; - } + psF32 slope; + + x.type.type = PS_TYPE_F32; + for (i=0;iimage->numRows;i++) { + for (j=0;jimage->numCols;j++) { + x.data.F32 = inputReadout->image->data.F32[i][j]; + binNum = p_psVectorBinDisect((psVector *)inFlux, &x); + + if (binNum == -2) { + // We get here if x is below the table lookup range. + inputReadout->image->data.F32[i][j] = outFlux->data.F32[0]; + numPixels++; + + } else if (binNum == -1) { + // We get here if x is above the table lookup range. + inputReadout->image->data.F32[i][j] = outFlux->data.F32[tableSize-1]; + numPixels++; + + } else if (binNum < -2) { + // We get here if there was some other problem. + psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): Could not perform p_psVectorBinDisect(). Returning inputReadout image."); + return(inputReadout); + numPixels++; + } else { + // Perform linear interpolation. + slope = (outFlux->data.F32[binNum+1] - outFlux->data.F32[binNum]) / + (inFlux->data.F32[binNum+1] - inFlux->data.F32[binNum]); + inputReadout->image->data.F32[i][j] = outFlux->data.F32[binNum] + + ((x.data.F32 - inFlux->data.F32[binNum]) * slope); } } - if (numPixels > 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels); - } } - + if (numPixels > 0) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels); + } return(inputReadout); } Index: trunk/psModules/src/detrend/pmNonLinear.h =================================================================== --- trunk/psModules/src/detrend/pmNonLinear.h (revision 5675) +++ trunk/psModules/src/detrend/pmNonLinear.h (revision 6872) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmNonLinear.h * @@ -5,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-20 23:06:24 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -16,15 +21,12 @@ #include "pslib.h" -#include "pmAstrometry.h" +#include "pmFPA.h" -pmReadout *pmNonLinearityPolynomial( - pmReadout *in, - const psPolynomial1D *coeff -); +pmReadout *pmNonLinearityPolynomial(pmReadout *in, + const psPolynomial1D *coeff); -pmReadout *pmNonLinearityLookup( - pmReadout *in, - const char *filename -); +pmReadout *pmNonLinearityLookup(pmReadout *in, + const psVector *inFlux, + const psVector *outFlux); #endif