Changeset 2585
- Timestamp:
- Dec 1, 2004, 11:26:17 AM (22 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/pmSubtractSky.c (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/pmSubtractSky.c
r2584 r2585 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1. 1$ $Name: not supported by cvs2svn $9 * @date $Date: 2004-12-01 21:2 0:24$8 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2004-12-01 21:26:17 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 32 32 } 33 33 34 PS_STAT_SAMPLE_MEAN = 0x000001, 35 PS_STAT_SAMPLE_MEDIAN = 0x000002, 36 PS_STAT_SAMPLE_STDEV = 0x000004, 37 PS_STAT_SAMPLE_QUARTILE = 0x000008, 38 PS_STAT_ROBUST_MEAN = 0x000010, 39 PS_STAT_ROBUST_MEDIAN = 0x000020, 40 PS_STAT_ROBUST_MODE = 0x000040, 41 PS_STAT_ROBUST_STDEV = 0x000080, 42 PS_STAT_ROBUST_QUARTILE = 0x000100, 43 PS_STAT_CLIPPED_MEAN = 0x000200, 44 PS_STAT_CLIPPED_STDEV = 0x000400, 45 PS_STAT_MAX = 0x000800, 46 PS_STAT_MIN = 0x001000, 47 PS_STAT_USE_RANGE = 0x002000, 48 PS_STAT_USE_BINSIZE = 0x004000, 49 PS_STAT_ROBUST_FOR_SAMPLE = 0x008000 50 51 psU64 getHighestPriorityStatOption(psU64 statOptions) 52 { 53 switch (statOptions) { 54 case PS_STAT_SAMPLE_MEAN: 55 return(PS_STAT_SAMPLE_MEAN); 56 case PS_STAT_SAMPLE_MEDIAN: 57 return(PS_STAT_SAMPLE_MEDIAN); 58 case PS_STAT_CLIPPED_MEAN: 59 return(PS_STAT_CLIPPED_MEAN); 60 case PS_STAT_ROBUST_MEAN: 61 return(PS_STAT_ROBUST_MEAN); 62 case PS_STAT_ROBUST_MEDIAN: 63 return(PS_STAT_ROBUST_MEDIAN); 64 case PS_STAT_ROBUST_MODE: 65 return(PS_STAT_ROBUST_MODE); 66 } 67 psError(XXX); 68 return(-1); 69 } 70 71 72 /****************************************************************************** 73 psImage *binImage(origImage, binFactor, statOptions): This routine takes an 74 input psImage and scales it smaller by a factor of binFactor. The statistic 75 used in combining input pixels is specified in statOptions. 76 77 XXX: use static vectors for myStats, binVector and binMask. 78 *****************************************************************************/ 79 psImage *binImage(psImage *origImage, 80 int binFactor, 81 psStatsOptions statOptions) 82 { 83 psVector *binVector = psVectorAlloc(binFactor * binFactor, PS_TYPE_F32); 84 psVector *binMask = psVectorAlloc(binFactor * binFactor, PS_TYPE_U8); 85 psStats *myStats = psStatsAlloc(statOptions) 86 87 if (binFactor <= 0) { 88 psLogMsg(__func__, PS_LOG_WARN, 89 "WARNING: binImage(): binFactor is %d\n", binFactor); 90 return(origImage); 91 } 92 if (binFactor == 1) { 93 return(origImage); 94 } 95 96 for (int row = 0; row < origImage->numRows ; row+=binFactor) { 97 for (int col = 0; col < origImage->numCols ; col+=binFactor) { 98 int count = 0; 99 for (int binRow = 0; binRow <= binFactor ; binRow++) { 100 for (int binCol = 0; binCol <= binFactor ; binCol++) { 101 if (((row + binRow) < origImage->numRows) && 102 ((col + binCol) < origImage->numCol)) { 103 binVector->data.F32[count] = 104 origImage->data.F32[row + binRow][col + binCol]; 105 binMask->data.U8[count] = 0; 106 } else { 107 binVector->data.F32[count] = 0.0; 108 binMask->data.U8[count] = 1; 109 } 110 count++; 111 } 112 } 113 myStats = psVectorStats(myStats, binVector, binMask, 1); 114 psF64 statValue; 115 psBool rc = p_psGetStatValue(myStats, &statValue); 116 117 if (rc == true) { 118 origImage->data.F32[row][col] = (psF32) statValue; 119 } else { 120 origImage->data.F32[row][col] = 0.0; 121 psLogMsg(__func__, PS_LOG_WARN, 122 "WARNING: pmSubtractSky(), binImage(): p_psGetStatValue() was FALSE\n"); 123 } 124 } 125 } 126 psFree(myStats); 127 psFree(binVector); 128 psFree(binMask); 129 130 return(origImage); 131 } 132 133 134 /****************************************************************************** 135 psReadout pmSubtractSky(): 136 137 XXX: use static vectors for myStats. 138 *****************************************************************************/ 139 psReadout *pmSubtractSky(psReadout *in, 140 void *fitSpec, 141 psFit fit, 142 int binFactor, 143 psStats *stats, 144 float clipSD) 145 { 146 if (fitSpec == NULL) { 147 return(in); 148 } 149 if ((fit != PM_FIT_POLYNOMIAL) && (fit != PM_FIT_SPLINE)) { 150 // No fit is specified. 151 return(in); 152 } 153 154 psImage *origImage = in->image; 155 psImage *binnedImage = NULL; 156 psPolynomial1D *myPoly; 157 psSpline1D *mySpline; 158 159 psStatsOptions statOptions = stats->options; 160 if (1 < p_psDetermineNumBits(statOptions) { 161 //XXX psWarning(PS_ERR_UNKNOWN,true, "Multiple statistical options have been requested.\n"); 162 statOptions = getHighestPriorityStatOption(statOptions) 163 ; 164 } 165 166 // Bin the input image according to input parameters. 167 if ((binFactor <= 0) || 168 (stats == NULL)) { 169 if (binFactor <= 0) { 170 psLogMsg(__func__, PS_LOG_WARN, 171 "WARNING: pmSubtractSky(): binFactor is %d\n", binFactor); 172 } 173 if (stats == NULL) { 174 psLogMsg(__func__, PS_LOG_WARN, 175 "WARNING: pmSubtractSky(): input parameter stats is NULL\n"); 176 } 177 binnedImage = origImage; 178 } else { 179 binnedImage = binImage(origImage, binFactor, statOptions); 180 } 181 182 // Clip pixels that are outside the acceptable range. 183 if (clipSD <= 0.0) { 184 psLogMsg(__func__, PS_LOG_WARN, 185 "WARNING: pmSubtractSky(): clipSD is %f\n", clipSD) 186 ; 187 } else { 188 psF64 binnedMean; 189 psF64 binnedStdev; 190 psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); 191 myStats = psImageStats(myStats, binedImage, NULL, 0); 192 p_psGetStatValue(myStats, &binnedMean); 193 194 myStats->options = PS_STAT_SAMPLE_STDEV; 195 myStats = psImageStats(myStats, binedImage, NULL, 0); 196 p_psGetStatValue(myStats, &binnedStdev); 197 psFree(myStats); 198 199 for (int row = 0; row < binnedImage->numRows ; row++) { 200 for (int col = 0; col < binnedImage->numCols ; col++) { 201 if (fabs(binnedImage->data.F32[row][col] - binnedMean) > 202 (clipSD * binnedStdev)) { 203 binnedImage->mask.U8[row][col] = 1; 204 } 205 } 206 } 207 } 208 209 // XXX: fit the polynomial to the binned image 210 if (fit == PM_FIT_POLYNOMIAL) { 211 // Fit a polynomial to the old overscan vector. 212 myPoly = (psPolynomial2D *) 213 fitSpec; 214 myPoly = psImageFitPolynomial(myPoly, origImage); 215 // XXX Do we need to do something with ordinate scaling if Chebyshev? 216 binnedImage = psImageEvalPolynomial(binnedImage, myPoly); 217 218 } else if (fit == PM_FIT_SPLINE) { 219 // Fit a spline to the old overscan vector. 220 mySpline = (psSpline1D *) 221 fitSpec; 222 // XXX: What do we do? We don't have 2-D splines. 223 return(in); 224 } 225 226 //Subtract the polynomially fitted image from the original image 227 if (binFactor <= 0) { 228 for (int row = 0; row < origImage->numRows ; row++) { 229 for (int col = 0; col < origImage->numCols ; col++) { 230 origImage->data.F32[row][col]-= binnedImage->data.F32[row][col]; 231 } 232 } 233 } else { 234 for (int row = 0; row < origImage->numRows ; row++) { 235 for (int col = 0; col < origImage->numCols ; col++) { 236 int binRow = row / binFactor; 237 int binCol = col / binFactor; 238 origImage->data.F32[row][col]-= binnedImage->data.F32[binRow][binCol]; 239 } 240 } 241 242 } 243 return(in); 244 } 34 35 36 psU64 getHighestPriorityStatOption(psU64 statOptions) 37 { 38 switch (statOptions) { 39 case PS_STAT_SAMPLE_MEAN: 40 return(PS_STAT_SAMPLE_MEAN); 41 case PS_STAT_SAMPLE_MEDIAN: 42 return(PS_STAT_SAMPLE_MEDIAN); 43 case PS_STAT_CLIPPED_MEAN: 44 return(PS_STAT_CLIPPED_MEAN); 45 case PS_STAT_ROBUST_MEAN: 46 return(PS_STAT_ROBUST_MEAN); 47 case PS_STAT_ROBUST_MEDIAN: 48 return(PS_STAT_ROBUST_MEDIAN); 49 case PS_STAT_ROBUST_MODE: 50 return(PS_STAT_ROBUST_MODE); 51 } 52 psError(XXX); 53 return(-1); 54 } 55 56 57 /****************************************************************************** 58 psImage *binImage(origImage, binFactor, statOptions): This routine takes an 59 input psImage and scales it smaller by a factor of binFactor. The statistic 60 used in combining input pixels is specified in statOptions. 61 62 XXX: use static vectors for myStats, binVector and binMask. 63 *****************************************************************************/ 64 psImage *binImage(psImage *origImage, 65 int binFactor, 66 psStatsOptions statOptions) 67 { 68 psVector *binVector = psVectorAlloc(binFactor * binFactor, PS_TYPE_F32); 69 psVector *binMask = psVectorAlloc(binFactor * binFactor, PS_TYPE_U8); 70 psStats *myStats = psStatsAlloc(statOptions) 71 72 if (binFactor <= 0) { 73 psLogMsg(__func__, PS_LOG_WARN, 74 "WARNING: binImage(): binFactor is %d\n", binFactor); 75 return(origImage); 76 } 77 if (binFactor == 1) { 78 return(origImage); 79 } 80 81 for (int row = 0; row < origImage->numRows ; row+=binFactor) { 82 for (int col = 0; col < origImage->numCols ; col+=binFactor) { 83 int count = 0; 84 for (int binRow = 0; binRow <= binFactor ; binRow++) { 85 for (int binCol = 0; binCol <= binFactor ; binCol++) { 86 if (((row + binRow) < origImage->numRows) && 87 ((col + binCol) < origImage->numCol)) { 88 binVector->data.F32[count] = 89 origImage->data.F32[row + binRow][col + binCol]; 90 binMask->data.U8[count] = 0; 91 } else { 92 binVector->data.F32[count] = 0.0; 93 binMask->data.U8[count] = 1; 94 } 95 count++; 96 } 97 } 98 myStats = psVectorStats(myStats, binVector, binMask, 1); 99 psF64 statValue; 100 psBool rc = p_psGetStatValue(myStats, &statValue); 101 102 if (rc == true) { 103 origImage->data.F32[row][col] = (psF32) statValue; 104 } else { 105 origImage->data.F32[row][col] = 0.0; 106 psLogMsg(__func__, PS_LOG_WARN, 107 "WARNING: pmSubtractSky(), binImage(): p_psGetStatValue() was FALSE\n"); 108 } 109 } 110 } 111 psFree(myStats); 112 psFree(binVector); 113 psFree(binMask); 114 115 return(origImage); 116 } 117 118 119 /****************************************************************************** 120 psReadout pmSubtractSky(): 121 122 XXX: use static vectors for myStats. 123 *****************************************************************************/ 124 psReadout *pmSubtractSky(psReadout *in, 125 void *fitSpec, 126 psFit fit, 127 int binFactor, 128 psStats *stats, 129 float clipSD) 130 { 131 if (fitSpec == NULL) { 132 return(in); 133 } 134 if ((fit != PM_FIT_POLYNOMIAL) && (fit != PM_FIT_SPLINE)) { 135 // No fit is specified. 136 return(in); 137 } 138 139 psImage *origImage = in->image; 140 psImage *binnedImage = NULL; 141 psPolynomial1D *myPoly; 142 psSpline1D *mySpline; 143 144 psStatsOptions statOptions = stats->options; 145 if (1 < p_psDetermineNumBits(statOptions) { 146 //XXX psWarning(PS_ERR_UNKNOWN,true, "Multiple statistical options have been requested.\n"); 147 statOptions = getHighestPriorityStatOption(statOptions) 148 ; 149 } 150 151 // Bin the input image according to input parameters. 152 if ((binFactor <= 0) || 153 (stats == NULL)) { 154 if (binFactor <= 0) { 155 psLogMsg(__func__, PS_LOG_WARN, 156 "WARNING: pmSubtractSky(): binFactor is %d\n", binFactor); 157 } 158 if (stats == NULL) { 159 psLogMsg(__func__, PS_LOG_WARN, 160 "WARNING: pmSubtractSky(): input parameter stats is NULL\n"); 161 } 162 binnedImage = origImage; 163 } else { 164 binnedImage = binImage(origImage, binFactor, statOptions); 165 } 166 167 // Clip pixels that are outside the acceptable range. 168 if (clipSD <= 0.0) { 169 psLogMsg(__func__, PS_LOG_WARN, 170 "WARNING: pmSubtractSky(): clipSD is %f\n", clipSD) 171 ; 172 } else { 173 psF64 binnedMean; 174 psF64 binnedStdev; 175 psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); 176 myStats = psImageStats(myStats, binedImage, NULL, 0); 177 p_psGetStatValue(myStats, &binnedMean); 178 179 myStats->options = PS_STAT_SAMPLE_STDEV; 180 myStats = psImageStats(myStats, binedImage, NULL, 0); 181 p_psGetStatValue(myStats, &binnedStdev); 182 psFree(myStats); 183 184 for (int row = 0; row < binnedImage->numRows ; row++) { 185 for (int col = 0; col < binnedImage->numCols ; col++) { 186 if (fabs(binnedImage->data.F32[row][col] - binnedMean) > 187 (clipSD * binnedStdev)) { 188 binnedImage->mask.U8[row][col] = 1; 189 } 190 } 191 } 192 } 193 194 // XXX: fit the polynomial to the binned image 195 if (fit == PM_FIT_POLYNOMIAL) { 196 // Fit a polynomial to the old overscan vector. 197 myPoly = (psPolynomial2D *) 198 fitSpec; 199 myPoly = psImageFitPolynomial(myPoly, origImage); 200 // XXX Do we need to do something with ordinate scaling if Chebyshev? 201 binnedImage = psImageEvalPolynomial(binnedImage, myPoly); 202 203 } else if (fit == PM_FIT_SPLINE) { 204 // Fit a spline to the old overscan vector. 205 mySpline = (psSpline1D *) 206 fitSpec; 207 // XXX: What do we do? We don't have 2-D splines. 208 return(in); 209 } 210 211 //Subtract the polynomially fitted image from the original image 212 if (binFactor <= 0) { 213 for (int row = 0; row < origImage->numRows ; row++) { 214 for (int col = 0; col < origImage->numCols ; col++) { 215 origImage->data.F32[row][col]-= binnedImage->data.F32[row][col]; 216 } 217 } 218 } else { 219 for (int row = 0; row < origImage->numRows ; row++) { 220 for (int col = 0; col < origImage->numCols ; col++) { 221 int binRow = row / binFactor; 222 int binCol = col / binFactor; 223 origImage->data.F32[row][col]-= binnedImage->data.F32[binRow][binCol]; 224 } 225 } 226 227 } 228 return(in); 229 }
Note:
See TracChangeset
for help on using the changeset viewer.
