Index: /branches/rel10_ifa/psModules/src/detrend/Makefile.am =================================================================== --- /branches/rel10_ifa/psModules/src/detrend/Makefile.am (revision 6841) +++ /branches/rel10_ifa/psModules/src/detrend/Makefile.am (revision 6842) @@ -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: /branches/rel10_ifa/psModules/src/detrend/pmFringeStats.c =================================================================== --- /branches/rel10_ifa/psModules/src/detrend/pmFringeStats.c (revision 6841) +++ /branches/rel10_ifa/psModules/src/detrend/pmFringeStats.c (revision 6842) @@ -3,6 +3,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-04 18:41:38 $ + * @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-12 23:50:45 $ * * Copyright 2004 IfA @@ -10,25 +10,22 @@ #include "pslib.h" +#include "pmFPA.h" #include "pmFringeStats.h" -static pmFringeStatsFree (pmFringeStats *fringe) -{ - psFree (fringe->x); - psFree (fringe->y); - psFree (fringe->f); - psFree (fringe->df); - psFree (fringe->mask); +static void fringeRegionsFree(pmFringeRegions *fringe) +{ + psFree(fringe->x); + psFree(fringe->y); + psFree(fringe->mask); return; } -pmFringeStats *pmFringeStatsAlloc (int nPts, int dX, int dY, int nX, int nY) -{ - pmFringeStats *fringe = psAlloc(sizeof(pmFringeStats)); - (void)psMemSetDeallocator(fringe, (psFreeFunc)pmFringeStatsFree); +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->f = NULL; - fringe->df = NULL; fringe->mask = NULL; @@ -44,24 +41,45 @@ } -static pmFringeScaleFree (pmFringeScale *scale) -{ - psFree (scale->coeff); - psFree (scale->coeffErr); + +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 *pmFringeScaleAlloc(int nFringeFrames) { pmFringeScale *scale = psAlloc(sizeof(pmFringeScale)); - (void)psMemSetDeallocator(scale, (psFreeFunc)pmFringeScaleFree); + (void)psMemSetDeallocator(scale, (psFreeFunc)fringeScaleFree); scale->nFringeFrames = nFringeFrames; - scale->coeff = psVectorAlloc (nFringeFrames + 1, PS_TYPE_F32); - scale->coeffErr = psVectorAlloc (nFringeFrames + 1, PS_TYPE_F32); + scale->coeff = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32); + scale->coeffErr = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32); return scale; } -bool pmFringeStatsCreatePoints (pmFringeStats *fringe, psImage *image) +bool pmFringeStatsCreatePoints(pmFringeRegions *fringe, psImage *image) { @@ -69,23 +87,24 @@ // 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); - - nX = image->numCols; - nY = image->numRows; + 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; - dX = fringe->dX; - dY = fringe->dY; + 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); + frnd = psRandomUniform(rnd); xPt[i] = (nX - 2*dX)* frnd + dX; - frnd = psRandomUniform (rnd); + frnd = psRandomUniform(rnd); yPt[i] = (nY - 2*dY)* frnd + dY; } @@ -93,26 +112,20 @@ } -bool pmFringeStatsMeasure(pmFringeStats *fringe, pmReadout *readout) -{ - if (fringe->x == NULL) { +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); - - psFree (fringe->f); - psFree (fringe->df); - psFree (fringe->mask); - - fringe->f = psVectorRecycle (fringe->f, fringe->x->n, PS_TYPE_F32); - fringe->df = psVectorRecycle (fringe->df, fringe->x->n, PS_TYPE_F32); - fringe->mask = psVectorRecycle (fringe->mask, fringe->x->n, PS_TYPE_U8); + 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 = fringe->f->data.F32; - psF32 *dfPt = fringe->df->data.F32; + psF32 *fPt = measurements->f->data.F32; + psF32 *dfPt = measurements->df->data.F32; int dX = fringe->dX; @@ -122,11 +135,12 @@ psImage *mask = readout->mask; - psStats *stats = psStatsAlloc (PS_STATS_MEDIAN | PS_STATS_STDEV); + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); for (int i = 0; i < fringe->x->n; i++) { - region = psRegionSet (xPt[i] - dX, yPt[i] - dY, xPt[i] + dX + 1, YPt[i] + dY + 1); - region = psRegionForImage (image, region); - raster = psImageSubset (image, region); - psImageStats (stats, image, mask, maskVal); + 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; @@ -134,34 +148,31 @@ } - return true; + return measurements; } // XXX include the fringe error (fringe->df) in the fit? -pmFringeScale *pmFringeScaleMeasure (pmFringeStats *science, psArray *fringes) -{ - +pmFringeScale *pmFringeScaleMeasure(pmFringeStats *science, psArray *fringes) +{ double sum; psVector *v1; psVector *v2; - pmFringeScale *scale = pmFringeScaleAlloc (fringes->n); + 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 (nCod, PS_TYPE_F64); + psImage *A = psImageAlloc(nCof, nCof, PS_TYPE_F64); + psVector *B = psVectorAlloc(nCof, PS_TYPE_F64); for (int i = 0; i < nCof; i++) { - if (i == 0) { - p1 = NULL; - } else { + psF32 *p1 = NULL; + if (i != 0) { v1 = fringes->data[i - 1]; p1 = v1->data.F32; } for (int j = 0; j < nCof; j++) { - if (j == 0) { - p2 = NULL; - } else { + psF32 *p2 = NULL; + if (j != 0) { v2 = fringes->data[j - 1]; p2 = v2->data.F32; @@ -169,7 +180,7 @@ sum = 0; - for (int I = 0; I < nPts; I++) { - f1 = (p1 == NULL) ? 1.0 : p1[I]; - f2 = (p2 == NULL) ? 1.0 : p2[I]; + 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; } @@ -178,8 +189,8 @@ sum = 0; - p2 = science->f->data.F32; - for (int I = 0; I < nPts; I++) { - f1 = (p1 == NULL) ? 1.0 : p1[I]; - f2 = p2[I]; + 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; } @@ -201,29 +212,29 @@ // XXX note that this modifies the input fringe images -psImage *pmFringeCorrect(psImage *out, psMetadata *info, psImage *science, psArray *fringeImage, psArray *fringeStats) -{ - - // measure the fringe stats - pmFringeScale *scale = pmFringeScaleMeasure (science, fringeStats); - - if (out != NULL) { - psImageInit (out, 0.0); - } +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 - for (int i = 0; i < fringeImage->n; i++) { + psImage *sumFringe = NULL; + for (int i = 0; i < fringeImages->n; i++) { // rescale the fringe image - psBinaryOp (fringeImage->data[i], fringeImage->data[i], "*", psScalarAlloc (scale->coeff->data.F32[i+1])); + psBinaryOp(fringeImages->data[i], fringeImages->data[i], "*", + psScalarAlloc(scale->coeff->data.F32[i+1], PS_TYPE_F32)); // sum together - out = psBinaryOp (out, out, "+", fringeImage->data[i]); - } + sumFringe = (psImage*)psBinaryOp(sumFringe, sumFringe, "+", fringeImages->data[i]); + } + psFree(scale); // subtract the resulting fringe frame - out = psBinaryOp (out, science, "-", out); - - return out; -} - + readout->image = (psImage*)psBinaryOp(readout->image, readout->image, "-", sumFringe); + + return sumFringe; +} + Index: /branches/rel10_ifa/psModules/src/detrend/pmFringeStats.h =================================================================== --- /branches/rel10_ifa/psModules/src/detrend/pmFringeStats.h (revision 6841) +++ /branches/rel10_ifa/psModules/src/detrend/pmFringeStats.h (revision 6842) @@ -5,6 +5,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-04 18:41:38 $ + * @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-12 23:50:45 $ * * Copyright 2004 IfA, University of Hawaii @@ -14,6 +14,5 @@ # define PM_FRINGE_STATS -/** Fringe Stats structure to hold the fringe statistics - * for a given image +/** Structure to hold the fringe measurement regions. */ typedef struct @@ -25,13 +24,11 @@ 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 f; // fringe point median - psVector df; // fringe point stdev - psVector mask; // fringe point on/off mask + psVector *x; // fringe point coordinates (col) + psVector *y; // fringe point coordinates (row) + psVector *mask; // fringe point on/off mask } -pmFringeStats; +pmFringeRegions; -pmFringeStats *pmFringeStatsAlloc ( +pmFringeRegions *pmFringeRegionsAlloc ( int nPts, // number of points to create int dX, // half-width of fringe boxes @@ -40,4 +37,17 @@ 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 @@ -60,7 +70,8 @@ * @return bool: True or false for success or failure */ -bool pmFringeStatsMeasure( - pmFringeStats *fringe, ///< results are placed here - psImage *image ///< measure fringes on this image +pmFringeStats *pmFringeStatsMeasure( + pmFringeRegions *fringe, ///< fringe regions + pmReadout *readout, ///< measure fringes on this readout + psMaskType maskVal ///< Mask value for statistics ); @@ -68,15 +79,14 @@ * * Given a science image, a set of fringe stats and a matched set of fringe images, correct the science image - * for the fringe images. + * for the fringe images. * - * @return psImage: corrected science image + * @return psImage: summed fringe image */ -psImage *pmFringeCorrect( - psImage *out, // place corrected image here (allocate if NULL) - psMetadata *info, // place info about results here - psImage *science, // image to be corrected - psArray *fringeImage, // fringe images to use in correction - psArray *fringeStats // fringe stats to use in correction -); +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