Changeset 5762
- Timestamp:
- Dec 12, 2005, 10:32:44 AM (21 years ago)
- Location:
- trunk/psModules/src/objects
- Files:
-
- 2 edited
-
pmPSFtry.c (modified) (5 diffs)
-
pmPSFtry.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/pmPSFtry.c
r5255 r5762 79 79 int Npsf = 0; 80 80 81 pmPSFtry *try 82 = pmPSFtryAlloc (sources, modelName); 81 pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName); 83 82 84 83 // stage 1: fit an independent model (freeModel) to all sources 85 84 psTimerStart ("fit"); 86 for (int i = 0; i < try 87 ->sources->n; i++) { 88 89 pmSource *source = try 90 ->sources->data[i]; 91 pmModel *model = pmSourceModelGuess (source, try 92 ->psf->type); 93 x = source->peak->x; 94 y = source->peak->y; 95 96 // set temporary object mask and fit object 97 // fit model as FLT, not PSF 98 psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); 99 status = pmSourceFitModel (source, model, false); 100 psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); 101 102 // exclude the poor fits 103 if (!status) { 104 try 105 ->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL; 106 psFree (model); 107 continue; 108 } 109 try 110 ->modelFLT->data[i] = model; 111 Nflt ++; 112 } 85 for (int i = 0; i < psfTry->sources->n; i++) { 86 87 pmSource *source = psfTry->sources->data[i]; 88 pmModel *model = pmSourceModelGuess (source, psfTry->psf->type); 89 x = source->peak->x; 90 y = source->peak->y; 91 92 // set temporary object mask and fit object 93 // fit model as FLT, not PSF 94 psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); 95 status = pmSourceFitModel (source, model, false); 96 psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); 97 98 // exclude the poor fits 99 if (!status) { 100 psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL; 101 psFree (model); 102 continue; 103 } 104 psfTry->modelFLT->data[i] = model; 105 Nflt ++; 106 } 113 107 psLogMsg ("psphot.psftry", 4, "fit flt: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); 114 108 psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n); 115 109 116 110 // make this optional? 117 // DumpModelFits ( try->modelFLT, "modelsFLT.dat");111 // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat"); 118 112 119 113 // stage 2: construct a psf (pmPSF) from this collection of model fits 120 pmPSFFromModels (try 121 ->psf, try 122 ->modelFLT, try 123 ->mask); 114 pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask); 124 115 125 116 // stage 3: refit with fixed shape parameters 126 117 psTimerStart ("fit"); 127 for (int i = 0; i < try 128 ->sources->n; i++) { 129 // masked for: bad model fit, outlier in parameters 130 if (try 131 ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; 132 133 pmSource *source = try 134 ->sources->data[i]; 135 pmModel *modelFLT = try 136 ->modelFLT->data[i]; 137 138 // set shape for this model based on PSF 139 pmModel *modelPSF = pmModelFromPSF (modelFLT, try 140 ->psf); 141 x = source->peak->x; 142 y = source->peak->y; 143 144 psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); 145 status = pmSourceFitModel (source, modelPSF, true); 146 147 // skip poor fits 148 if (!status) { 149 try 150 ->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL; 151 psFree (modelPSF); 152 goto next_source; 153 } 154 155 // otherwise, save the resulting model 156 try 157 ->modelPSF->data[i] = modelPSF; 158 159 // XXX : use a different aperture radius from the fit radius? 160 // XXX : use a different estimator for the local sky? 161 // XXX : pass 'source' as input? 162 if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) { 163 try 164 ->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; 165 goto next_source; 166 } 167 168 try 169 ->metric->data.F64[i] = obsMag - fitMag; 170 try 171 ->fitMag->data.F64[i] = fitMag; 172 Npsf ++; 118 for (int i = 0; i < psfTry->sources->n; i++) { 119 // masked for: bad model fit, outlier in parameters 120 if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL) 121 continue; 122 123 pmSource *source = psfTry->sources->data[i]; 124 pmModel *modelFLT = psfTry->modelFLT->data[i]; 125 126 // set shape for this model based on PSF 127 pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf); 128 x = source->peak->x; 129 y = source->peak->y; 130 131 psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); 132 status = pmSourceFitModel (source, modelPSF, true); 133 134 // skip poor fits 135 if (!status) { 136 psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL; 137 psFree (modelPSF); 138 goto next_source; 139 } 140 141 // otherwise, save the resulting model 142 psfTry->modelPSF->data[i] = modelPSF; 143 144 // XXX : use a different aperture radius from the fit radius? 145 // XXX : use a different estimator for the local sky? 146 // XXX : pass 'source' as input? 147 if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) { 148 psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; 149 goto next_source; 150 } 151 152 psfTry->metric->data.F64[i] = obsMag - fitMag; 153 psfTry->fitMag->data.F64[i] = fitMag; 154 Npsf ++; 173 155 174 156 next_source: 175 psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);176 177 }157 psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); 158 159 } 178 160 psLogMsg ("psphot.psftry", 4, "fit psf: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); 179 161 psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n); 180 162 181 163 // make this optional 182 // DumpModelFits ( try->modelPSF, "modelsPSF.dat");164 // DumpModelFits (psfTry->modelPSF, "modelsPSF.dat"); 183 165 184 166 // XXX this function wants aperture radius for pmSourcePhotometry 185 pmPSFtryMetric (try 186 , RADIUS); 167 pmPSFtryMetric (psfTry, RADIUS); 187 168 psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n", 188 modelName, try 189 ->psf->ApResid, try 190 ->psf->dApResid, try 191 ->psf->skyBias); 192 193 return (try 194 ); 169 modelName, 170 psfTry->psf->ApResid, 171 psfTry->psf->dApResid, 172 psfTry->psf->skyBias); 173 174 return (psfTry); 195 175 } 196 176 197 177 198 bool pmPSFtryMetric (pmPSFtry *try 199 , float RADIUS) 178 bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS) 200 179 { 201 180 … … 203 182 int nKeep, nSkip; 204 183 205 // the measured (aperture - fit) magnitudes (dA == try->metric)184 // the measured (aperture - fit) magnitudes (dA == psfTry->metric) 206 185 // depend on both the true ap-fit (dAo) and the bias in the sky measurement: 207 186 // dA = dAo + dsky/flux … … 213 192 214 193 // rflux = ten(0.4*fitMag); 215 psVector *rflux = psVectorAlloc (try 216 ->sources->n, PS_TYPE_F64); 217 for (int i = 0; i < try 218 ->sources->n; i++) { 219 if (try 220 ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; 221 rflux->data.F64[i] = pow(10.0, 0.4*try 222 ->fitMag->data.F64[i]); 223 } 194 psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); 195 for (int i = 0; i < psfTry->sources->n; i++) { 196 if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL) 197 continue; 198 rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]); 199 } 224 200 225 201 // find min and max of (1/flux): 226 202 psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); 227 psVectorStats (stats, rflux, NULL, try 228 ->mask, PSFTRY_MASK_ALL); 203 psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL); 229 204 230 205 // build binned versions of rflux, metric … … 240 215 for (int i = 0; i < daBin->n; i++) { 241 216 242 psVector *tmp = psVectorAlloc (try 243 ->sources->n, PS_TYPE_F64); 217 psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); 244 218 tmp->n = 0; 245 219 246 220 // accumulate data within bin range 247 for (int j = 0; j < try 248 ->sources->n; j++) { 249 // masked for: bad model fit, outlier in parameters 250 if (try 251 ->mask->data.U8[j] & PSFTRY_MASK_ALL) continue; 252 253 // skip points with extreme dA values 254 if (fabs(try 255 ->metric->data.F64[j]) > 0.5) continue; 256 257 // skip points outside of this bin 258 if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) 259 continue; 260 if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) 261 continue; 262 263 tmp->data.F64[tmp->n] = try 264 ->metric->data.F64[j]; 265 tmp->n ++; 266 } 221 for (int j = 0; j < psfTry->sources->n; j++) { 222 // masked for: bad model fit, outlier in parameters 223 if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL) 224 continue; 225 226 // skip points with extreme dA values 227 if (fabs(psfTry->metric->data.F64[j]) > 0.5) 228 continue; 229 230 // skip points outside of this bin 231 if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) 232 continue; 233 if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) 234 continue; 235 236 tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j]; 237 tmp->n ++; 238 } 267 239 268 240 // is this a valid point? … … 319 291 stats = psVectorStats (stats, daResid, NULL, maskB, 1); 320 292 321 try 322 ->psf->ApResid = poly->coeff[0]; 323 try 324 ->psf->dApResid = stats->clippedStdev; 325 try 326 ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); 293 psfTry->psf->ApResid = poly->coeff[0]; 294 psfTry->psf->dApResid = stats->clippedStdev; 295 psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); 327 296 328 297 psFree (rflux); -
trunk/psModules/src/objects/pmPSFtry.h
r5255 r5762 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1. 1$ $Name: not supported by cvs2svn $9 * @date $Date: 2005-1 0-10 19:53:40$8 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2005-12-12 20:32:44 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 105 105 */ 106 106 bool pmPSFtryMetric( 107 pmPSFtry *try 108 , ///< Add comment. 109 float RADIUS ///< Add comment. 107 pmPSFtry *psfTry, ///< Add comment. 108 float RADIUS ///< Add comment. 110 109 ); 111 110
Note:
See TracChangeset
for help on using the changeset viewer.
