Changeset 6842
- Timestamp:
- Apr 12, 2006, 1:50:45 PM (20 years ago)
- Location:
- branches/rel10_ifa/psModules/src/detrend
- Files:
-
- 3 edited
-
Makefile.am (modified) (1 diff)
-
pmFringeStats.c (modified) (10 diffs)
-
pmFringeStats.h (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/rel10_ifa/psModules/src/detrend/Makefile.am
r5170 r6842 3 3 libpsmoduledetrend_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS) 4 4 libpsmoduledetrend_la_LDFLAGS = -release $(PACKAGE_VERSION) 5 libpsmoduledetrend_la_SOURCES = pmFlatField.c \ 6 pmMaskBadPixels.c \ 7 pmNonLinear.c 5 libpsmoduledetrend_la_SOURCES = \ 6 pmFlatField.c \ 7 pmFringeStats.c \ 8 pmMaskBadPixels.c \ 9 pmNonLinear.c 8 10 9 11 psmoduleincludedir = $(includedir) 10 12 psmoduleinclude_HEADERS = \ 11 pmFlatField.h \ 12 pmFlatFieldErrors.h \ 13 pmMaskBadPixelsErrors.h \ 14 pmMaskBadPixels.h \ 15 pmNonLinear.h 13 pmFlatField.h \ 14 pmFlatFieldErrors.h \ 15 pmFringeStats.h \ 16 pmMaskBadPixelsErrors.h \ 17 pmMaskBadPixels.h \ 18 pmNonLinear.h 16 19 17 20 EXTRA_DIST = pmFlatFieldErrors.dat pmMaskBadPixelsErrors.dat -
branches/rel10_ifa/psModules/src/detrend/pmFringeStats.c
r6761 r6842 3 3 * @author Eugene Magnier, IfA 4 4 * 5 * @version $Revision: 1.1.2. 2$ $Name: not supported by cvs2svn $6 * @date $Date: 2006-04- 04 18:41:38$5 * @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2006-04-12 23:50:45 $ 7 7 * 8 8 * Copyright 2004 IfA … … 10 10 11 11 #include "pslib.h" 12 #include "pmFPA.h" 12 13 #include "pmFringeStats.h" 13 14 14 static pmFringeStatsFree (pmFringeStats *fringe) 15 { 16 psFree (fringe->x); 17 psFree (fringe->y); 18 psFree (fringe->f); 19 psFree (fringe->df); 20 psFree (fringe->mask); 15 static void fringeRegionsFree(pmFringeRegions *fringe) 16 { 17 psFree(fringe->x); 18 psFree(fringe->y); 19 psFree(fringe->mask); 21 20 return; 22 21 } 23 22 24 pmFringe Stats *pmFringeStatsAlloc(int nPts, int dX, int dY, int nX, int nY)25 { 26 pmFringe Stats *fringe = psAlloc(sizeof(pmFringeStats));27 (void)psMemSetDeallocator(fringe, (psFreeFunc) pmFringeStatsFree);23 pmFringeRegions *pmFringeRegionsAlloc(int nPts, int dX, int dY, int nX, int nY) 24 { 25 pmFringeRegions *fringe = psAlloc(sizeof(pmFringeRegions)); 26 (void)psMemSetDeallocator(fringe, (psFreeFunc)fringeRegionsFree); 28 27 29 28 fringe->x = NULL; 30 29 fringe->y = NULL; 31 fringe->f = NULL;32 fringe->df = NULL;33 30 fringe->mask = NULL; 34 31 … … 44 41 } 45 42 46 static pmFringeScaleFree (pmFringeScale *scale) 47 { 48 psFree (scale->coeff); 49 psFree (scale->coeffErr); 43 44 static void fringeStatsFree(pmFringeStats *stats) 45 { 46 psFree(stats->regions); 47 psFree(stats->f); 48 psFree(stats->df); 49 } 50 51 pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions) 52 { 53 pmFringeStats *stats = psAlloc(sizeof(pmFringeStats)); 54 (void)psMemSetDeallocator(stats, (psFreeFunc)fringeStatsFree); 55 56 stats->regions = psMemIncrRefCounter(regions); 57 stats->f = psVectorAlloc(regions->x->n, PS_TYPE_F32); 58 stats->df = psVectorAlloc(regions->x->n, PS_TYPE_F32); 59 60 return stats; 61 } 62 63 64 static void fringeScaleFree(pmFringeScale *scale) 65 { 66 psFree(scale->coeff); 67 psFree(scale->coeffErr); 50 68 return; 51 69 } 52 70 53 pmFringeScale *pmFringeScaleAlloc (int nFringeFrames)71 pmFringeScale *pmFringeScaleAlloc(int nFringeFrames) 54 72 { 55 73 pmFringeScale *scale = psAlloc(sizeof(pmFringeScale)); 56 (void)psMemSetDeallocator(scale, (psFreeFunc) pmFringeScaleFree);74 (void)psMemSetDeallocator(scale, (psFreeFunc)fringeScaleFree); 57 75 58 76 scale->nFringeFrames = nFringeFrames; 59 scale->coeff = psVectorAlloc (nFringeFrames + 1, PS_TYPE_F32);60 scale->coeffErr = psVectorAlloc (nFringeFrames + 1, PS_TYPE_F32);77 scale->coeff = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32); 78 scale->coeffErr = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32); 61 79 62 80 return scale; 63 81 } 64 82 65 bool pmFringeStatsCreatePoints (pmFringeStats *fringe, psImage *image)83 bool pmFringeStatsCreatePoints(pmFringeRegions *fringe, psImage *image) 66 84 { 67 85 … … 69 87 // create fringe->nRequested 70 88 71 psRandom *rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); 72 73 fringe->x = psVectorRecycle (fringe->x, fringe->nRequested, PS_TYPE_F32); 74 fringe->y = psVectorRecycle (fringe->y, fringe->nRequested, PS_TYPE_F32); 75 76 nX = image->numCols; 77 nY = image->numRows; 89 psRandom *rnd = psRandomAlloc(PS_RANDOM_TAUS, 0); 90 91 fringe->x = psVectorRecycle(fringe->x, fringe->nRequested, PS_TYPE_F32); 92 fringe->y = psVectorRecycle(fringe->y, fringe->nRequested, PS_TYPE_F32); 93 fringe->mask = psVectorRecycle(fringe->mask, fringe->nRequested, PS_TYPE_U8); 94 95 int nX = image->numCols; 96 int nY = image->numRows; 78 97 79 98 psF32 *xPt = fringe->x->data.F32; 80 99 psF32 *yPt = fringe->y->data.F32; 81 100 82 dX = fringe->dX;83 dY = fringe->dY;101 int dX = fringe->dX; 102 int dY = fringe->dY; 84 103 85 104 // generate random points located within image bounds 86 105 for (int i = 0; i < fringe->nRequested; i++) { 87 frnd = psRandomUniform (rnd);106 frnd = psRandomUniform(rnd); 88 107 xPt[i] = (nX - 2*dX)* frnd + dX; 89 frnd = psRandomUniform (rnd);108 frnd = psRandomUniform(rnd); 90 109 yPt[i] = (nY - 2*dY)* frnd + dY; 91 110 } … … 93 112 } 94 113 95 bool pmFringeStatsMeasure(pmFringeStats *fringe, pmReadout *readout)96 { 97 if ( fringe->x == NULL) {114 pmFringeStats *pmFringeStatsMeasure(pmFringeRegions *fringe, pmReadout *readout, psMaskType maskVal) 115 { 116 if (!fringe->x || !fringe->y) { 98 117 // create the fringe vectors for this image 99 pmFringeStatsCreatePoints (fringe, readout->image); 100 } 101 102 PS_ASSERT_PTR_NON_NULL (fringe->x, false); 103 PS_ASSERT_PTR_NON_NULL (fringe->y, false); 104 105 psFree (fringe->f); 106 psFree (fringe->df); 107 psFree (fringe->mask); 108 109 fringe->f = psVectorRecycle (fringe->f, fringe->x->n, PS_TYPE_F32); 110 fringe->df = psVectorRecycle (fringe->df, fringe->x->n, PS_TYPE_F32); 111 fringe->mask = psVectorRecycle (fringe->mask, fringe->x->n, PS_TYPE_U8); 118 pmFringeStatsCreatePoints(fringe, readout->image); 119 } 120 121 PS_ASSERT_PTR_NON_NULL(fringe->x, false); 122 PS_ASSERT_PTR_NON_NULL(fringe->y, false); 123 124 pmFringeStats *measurements = pmFringeStatsAlloc(fringe); 112 125 113 126 psF32 *xPt = fringe->x->data.F32; 114 127 psF32 *yPt = fringe->y->data.F32; 115 psF32 *fPt = fringe->f->data.F32;116 psF32 *dfPt = fringe->df->data.F32;128 psF32 *fPt = measurements->f->data.F32; 129 psF32 *dfPt = measurements->df->data.F32; 117 130 118 131 int dX = fringe->dX; … … 122 135 psImage *mask = readout->mask; 123 136 124 psStats *stats = psStatsAlloc (PS_STATS_MEDIAN | PS_STATS_STDEV);137 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); 125 138 126 139 for (int i = 0; i < fringe->x->n; i++) { 127 region = psRegionSet (xPt[i] - dX, yPt[i] - dY, xPt[i] + dX + 1, YPt[i] + dY + 1); 128 region = psRegionForImage (image, region); 129 raster = psImageSubset (image, region); 130 psImageStats (stats, image, mask, maskVal); 140 psRegion region = psRegionSet(xPt[i] - dX, yPt[i] - dY, xPt[i] + dX + 1, yPt[i] + dY + 1); 141 region = psRegionForImage(image, region); 142 psImage *subImage = psImageSubset(image, region); 143 psImage *subMask = psImageSubset(mask, region); 144 psImageStats(stats, subImage, subMask, maskVal); 131 145 132 146 fPt[i] = stats->sampleMedian; … … 134 148 } 135 149 136 return true;150 return measurements; 137 151 } 138 152 139 153 // XXX include the fringe error (fringe->df) in the fit? 140 pmFringeScale *pmFringeScaleMeasure (pmFringeStats *science, psArray *fringes) 141 { 142 154 pmFringeScale *pmFringeScaleMeasure(pmFringeStats *science, psArray *fringes) 155 { 143 156 double sum; 144 157 psVector *v1; 145 158 psVector *v2; 146 159 147 pmFringeScale *scale = pmFringeScaleAlloc (fringes->n);160 pmFringeScale *scale = pmFringeScaleAlloc(fringes->n); 148 161 149 162 int nCof = fringes->n + 1; 150 163 int nPts = science->f->n; 151 164 152 psImage *A = psImageAlloc (nCof, nCof, PS_TYPE_F64);153 psVector *B = psVectorAlloc (nCod, PS_TYPE_F64);165 psImage *A = psImageAlloc(nCof, nCof, PS_TYPE_F64); 166 psVector *B = psVectorAlloc(nCof, PS_TYPE_F64); 154 167 155 168 for (int i = 0; i < nCof; i++) { 156 if (i == 0) { 157 p1 = NULL; 158 } else { 169 psF32 *p1 = NULL; 170 if (i != 0) { 159 171 v1 = fringes->data[i - 1]; 160 172 p1 = v1->data.F32; 161 173 } 162 174 for (int j = 0; j < nCof; j++) { 163 if (j == 0) { 164 p2 = NULL; 165 } else { 175 psF32 *p2 = NULL; 176 if (j != 0) { 166 177 v2 = fringes->data[j - 1]; 167 178 p2 = v2->data.F32; … … 169 180 170 181 sum = 0; 171 for (int I = 0; I < nPts; I++) {172 f1 = (p1 == NULL) ? 1.0 : p1[I];173 f2 = (p2 == NULL) ? 1.0 : p2[I];182 for (int k = 0; k < nPts; k++) { 183 psF32 f1 = (p1 == NULL) ? 1.0 : p1[k]; 184 psF32 f2 = (p2 == NULL) ? 1.0 : p2[k]; 174 185 sum += f1*f2; 175 186 } … … 178 189 179 190 sum = 0; 180 p 2 = science->f->data.F32;181 for (int I = 0; I < nPts; I++) {182 f1 = (p1 == NULL) ? 1.0 : p1[I];183 f2 = p2[I];191 psF32 *p2 = science->f->data.F32; 192 for (int j = 0; j < nPts; j++) { 193 psF32 f1 = (p1 == NULL) ? 1.0 : p1[j]; 194 psF32 f2 = p2[j]; 184 195 sum += f1*f2; 185 196 } … … 201 212 202 213 // XXX note that this modifies the input fringe images 203 psImage *pmFringeCorrect(psImage *out, psMetadata *info, psImage *science, psArray *fringeImage, psArray *fringeStats) 204 { 205 206 // measure the fringe stats 207 pmFringeScale *scale = pmFringeScaleMeasure (science, fringeStats); 208 209 if (out != NULL) { 210 psImageInit (out, 0.0); 211 } 214 psImage *pmFringeCorrect(pmReadout *readout, pmFringeRegions *fringes, psArray *fringeImages, psArray *fringeStats, psMaskType maskVal) 215 { 216 // measure the fringe stats for the science frame and solve for the scales 217 pmFringeStats *scienceStats = pmFringeStatsMeasure(fringes, readout, maskVal); 218 pmFringeScale *scale = pmFringeScaleMeasure(scienceStats, fringeStats); 219 psFree(scienceStats); 212 220 213 221 // build the fringe correction image 214 222 // XXX we could save data space by making the first image the output image 215 for (int i = 0; i < fringeImage->n; i++) { 223 psImage *sumFringe = NULL; 224 for (int i = 0; i < fringeImages->n; i++) { 216 225 217 226 // rescale the fringe image 218 psBinaryOp (fringeImage->data[i], fringeImage->data[i], "*", psScalarAlloc (scale->coeff->data.F32[i+1])); 227 psBinaryOp(fringeImages->data[i], fringeImages->data[i], "*", 228 psScalarAlloc(scale->coeff->data.F32[i+1], PS_TYPE_F32)); 219 229 220 230 // sum together 221 out = psBinaryOp (out, out, "+", fringeImage->data[i]); 222 } 231 sumFringe = (psImage*)psBinaryOp(sumFringe, sumFringe, "+", fringeImages->data[i]); 232 } 233 psFree(scale); 223 234 224 235 // subtract the resulting fringe frame 225 out = psBinaryOp (out, science, "-", out);226 227 return out;228 } 229 236 readout->image = (psImage*)psBinaryOp(readout->image, readout->image, "-", sumFringe); 237 238 return sumFringe; 239 } 240 -
branches/rel10_ifa/psModules/src/detrend/pmFringeStats.h
r6761 r6842 5 5 * @author Eugene Magnier, IfA 6 6 * 7 * @version $Revision: 1.1.2. 2$ $Name: not supported by cvs2svn $8 * @date $Date: 2006-04- 04 18:41:38$7 * @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2006-04-12 23:50:45 $ 9 9 * 10 10 * Copyright 2004 IfA, University of Hawaii … … 14 14 # define PM_FRINGE_STATS 15 15 16 /** Fringe Stats structure to hold the fringe statistics 17 * for a given image 16 /** Structure to hold the fringe measurement regions. 18 17 */ 19 18 typedef struct … … 25 24 psU32 nX; // large-scale smoothing in x (col) 26 25 psU32 nY; // large-scale smoothing in y (row) 27 psVector x; // fringe point coordinates (col) 28 psVector y; // fringe point coordinates (row) 29 psVector f; // fringe point median 30 psVector df; // fringe point stdev 31 psVector mask; // fringe point on/off mask 26 psVector *x; // fringe point coordinates (col) 27 psVector *y; // fringe point coordinates (row) 28 psVector *mask; // fringe point on/off mask 32 29 } 33 pmFringe Stats;30 pmFringeRegions; 34 31 35 pmFringe Stats *pmFringeStatsAlloc (32 pmFringeRegions *pmFringeRegionsAlloc ( 36 33 int nPts, // number of points to create 37 34 int dX, // half-width of fringe boxes … … 40 37 int nY // smoothing scale in y 41 38 ); 39 40 /** Structure to hold the fringe measurements for a particular image 41 */ 42 typedef struct 43 { 44 pmFringeRegions *regions; // Fringe regions 45 psVector *f; // fringe point median 46 psVector *df; // fringe point stdev 47 } 48 pmFringeStats; 49 50 pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions); 51 42 52 43 53 /** the pmFringeScale structure defines the relationship between two fringe measurements … … 60 70 * @return bool: True or false for success or failure 61 71 */ 62 bool pmFringeStatsMeasure( 63 pmFringeStats *fringe, ///< results are placed here 64 psImage *image ///< measure fringes on this image 72 pmFringeStats *pmFringeStatsMeasure( 73 pmFringeRegions *fringe, ///< fringe regions 74 pmReadout *readout, ///< measure fringes on this readout 75 psMaskType maskVal ///< Mask value for statistics 65 76 ); 66 77 … … 68 79 * 69 80 * Given a science image, a set of fringe stats and a matched set of fringe images, correct the science image 70 * for the fringe images. 81 * for the fringe images. 71 82 * 72 * @return psImage: corrected science image83 * @return psImage: summed fringe image 73 84 */ 74 psImage *pmFringeCorrect( 75 psImage *out, // place corrected image here (allocate if NULL) 76 psMetadata *info, // place info about results here 77 psImage *science, // image to be corrected 78 psArray *fringeImage, // fringe images to use in correction 79 psArray *fringeStats // fringe stats to use in correction 80 ); 85 psImage *pmFringeCorrect(pmReadout *in, // place info about results here 86 pmFringeRegions *fringes, // The fringe regions used 87 psArray *fringeImages, // fringe images to use in correction 88 psArray *fringeStats, // fringe stats to use in correction 89 psMaskType maskVal // Value to mask 90 ); 81 91 82 92 # endif
Note:
See TracChangeset
for help on using the changeset viewer.
