Changeset 9734 for trunk/psphot/src/psphotApResid.c
- Timestamp:
- Oct 24, 2006, 2:07:03 PM (20 years ago)
- File:
-
- 1 edited
-
trunk/psphot/src/psphotApResid.c (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot/src/psphotApResid.c
r9568 r9734 4 4 5 5 // measure the aperture residual statistics 6 bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { 6 bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { 7 7 8 8 int Nfail = 0; … … 18 18 float MAX_AP_OFFSET = psMetadataLookupF32 (&status, recipe, "MAX_AP_OFFSET"); 19 19 20 // set limits on the aperture magnitudes 20 // set limits on the aperture magnitudes 21 21 pmSourceMagnitudesInit (recipe); 22 22 … … 25 25 psf->growth = pmGrowthCurveAlloc (3.0, 100.0, REF_RADIUS); 26 26 psphotGrowthCurve (readout, psf); 27 28 psVector *mask = psVectorAlloc (300, PS_TYPE_U8); 29 psVector *xPos = psVectorAlloc (300, PS_TYPE_F64); 30 psVector *yPos = psVectorAlloc (300, PS_TYPE_F64); 31 psVector *flux = psVectorAlloc (300, PS_TYPE_F64); 32 psVector *r2rflux = psVectorAlloc (300, PS_TYPE_F64); 33 psVector *apResid = psVectorAlloc (300, PS_TYPE_F64); 34 psVector *dMag = psVectorAlloc (300, PS_TYPE_F64); 35 mask->n = xPos->n = yPos->n = flux->n = r2rflux->n = apResid->n = dMag->n = 0; 27 28 psVector *mask = psVectorAllocEmpty (300, PS_TYPE_U8); 29 psVector *xPos = psVectorAllocEmpty (300, PS_TYPE_F64); 30 psVector *yPos = psVectorAllocEmpty (300, PS_TYPE_F64); 31 psVector *flux = psVectorAllocEmpty (300, PS_TYPE_F64); 32 psVector *r2rflux = psVectorAllocEmpty (300, PS_TYPE_F64); 33 psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F64); 34 psVector *dMag = psVectorAllocEmpty (300, PS_TYPE_F64); 36 35 Npsf = 0; 37 36 38 37 // select all good PM_SOURCE_TYPE_STAR entries 39 38 for (int i = 0; i < sources->n; i++) { 40 source = sources->data[i];41 model = source->modelPSF;42 43 if (source->type != PM_SOURCE_TYPE_STAR) continue; 44 if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;45 if (source->mode & PM_SOURCE_MODE_BLEND) continue;46 if (source->mode & PM_SOURCE_MODE_FAIL) continue;47 if (source->mode & PM_SOURCE_MODE_POOR) continue;48 49 // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures50 // will fail if below S/N threshold or model is missing51 if (!pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH)) {52 Nskip ++;53 continue;54 }55 56 apResid->data.F64[Npsf] = source->apMag - source->psfMag;57 58 xPos->data.F64[Npsf] = model->params->data.F32[PM_PAR_XPOS];59 yPos->data.F64[Npsf] = model->params->data.F32[PM_PAR_YPOS];60 61 flux->data.F64[Npsf] = pow(10.0, -0.4*source->psfMag);62 r2rflux->data.F64[Npsf] = PS_SQR(model->radiusTMP) / flux->data.F64[Npsf];63 64 mask->data.U8[Npsf] = 0;65 66 // XXX sanity clip?67 // XXX need to see if all data were tossed?68 // XXX need to subtract median?69 // XXX need to put this in the config data...70 // if (fabs(apResid->data.F64[Npsf]) > 0.2) continue;71 if ((MAX_AP_OFFSET > 0) && (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET)) {72 Nfail ++;73 continue;74 }75 76 dMag->data.F64[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];77 78 psVectorExtend (mask,100, 1);79 psVectorExtend (xPos,100, 1);80 psVectorExtend (yPos,100, 1);81 psVectorExtend (flux, 100, 1);82 psVectorExtend (r2rflux, 100, 1);83 psVectorExtend (dMag, 100, 1);84 psVectorExtend (apResid, 100, 1);85 Npsf ++;86 } 87 psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects (%d skipped, %d failed, %ld invalid)\n", 88 psTimerMark ("psphot"), Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);39 source = sources->data[i]; 40 model = source->modelPSF; 41 42 if (source->type != PM_SOURCE_TYPE_STAR) continue; 43 if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; 44 if (source->mode & PM_SOURCE_MODE_BLEND) continue; 45 if (source->mode & PM_SOURCE_MODE_FAIL) continue; 46 if (source->mode & PM_SOURCE_MODE_POOR) continue; 47 48 // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures 49 // will fail if below S/N threshold or model is missing 50 if (!pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH)) { 51 Nskip ++; 52 continue; 53 } 54 55 apResid->data.F64[Npsf] = source->apMag - source->psfMag; 56 57 xPos->data.F64[Npsf] = model->params->data.F32[PM_PAR_XPOS]; 58 yPos->data.F64[Npsf] = model->params->data.F32[PM_PAR_YPOS]; 59 60 flux->data.F64[Npsf] = pow(10.0, -0.4*source->psfMag); 61 r2rflux->data.F64[Npsf] = PS_SQR(model->radiusTMP) / flux->data.F64[Npsf]; 62 63 mask->data.U8[Npsf] = 0; 64 65 // XXX sanity clip? 66 // XXX need to see if all data were tossed? 67 // XXX need to subtract median? 68 // XXX need to put this in the config data... 69 // if (fabs(apResid->data.F64[Npsf]) > 0.2) continue; 70 if ((MAX_AP_OFFSET > 0) && (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET)) { 71 Nfail ++; 72 continue; 73 } 74 75 dMag->data.F64[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; 76 77 psVectorExtend (mask, 100, 1); 78 psVectorExtend (xPos, 100, 1); 79 psVectorExtend (yPos, 100, 1); 80 psVectorExtend (flux, 100, 1); 81 psVectorExtend (r2rflux, 100, 1); 82 psVectorExtend (dMag, 100, 1); 83 psVectorExtend (apResid, 100, 1); 84 Npsf ++; 85 } 86 psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects (%d skipped, %d failed, %ld invalid)\n", 87 psTimerMark ("psphot"), Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail); 89 88 90 89 // XXX choose a better value here? 91 90 if (Npsf < 5) { 92 psError(PSPHOT_ERR_APERTURE, true, "no valid aperture residual sources, giving up");93 return false;94 } 95 96 // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL 91 psError(PSPHOT_ERR_APERTURE, true, "no valid aperture residual sources, giving up"); 92 return false; 93 } 94 95 // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL 97 96 // APTREND options are used in the switch block below 98 97 pmPSFApTrendOptions ApTrendOption = DEFAULT_OPTION; … … 100 99 if (status) ApTrendOption = pmPSFApTrendOptionFromName (optionName); 101 100 if (ApTrendOption == PM_PSF_APTREND_ERROR) { 102 psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName);103 return false;101 psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName); 102 return false; 104 103 } 105 104 … … 112 111 switch (ApTrendOption) { 113 112 case PM_PSF_APTREND_NONE: 114 // remove ApTrend fit from pmPSFtry115 psf->ApTrend->coeff[0][0][0][0] = 0;116 break;113 // remove ApTrend fit from pmPSFtry 114 psf->ApTrend->coeff[0][0][0][0] = 0; 115 break; 117 116 case PM_PSF_APTREND_CONSTANT: 118 stats->clipIter = 2;119 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);120 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);121 if (psf->ApTrend == NULL) {122 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");123 return false;124 }125 // apply the fit126 stats->clipIter = 3;127 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);128 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);129 if (psf->ApTrend == NULL) {130 psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");131 return false;132 }133 break;117 stats->clipIter = 2; 118 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); 119 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 120 if (psf->ApTrend == NULL) { 121 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); 122 return false; 123 } 124 // apply the fit 125 stats->clipIter = 3; 126 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); 127 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 128 if (psf->ApTrend == NULL) { 129 psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction"); 130 return false; 131 } 132 break; 134 133 case PM_PSF_APTREND_SKYBIAS: 135 // first clip out objects which are too far from the median 136 stats->clipIter = 2;137 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);138 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);139 if (psf->ApTrend == NULL) {140 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");141 return false;142 }143 // apply the fit144 stats->clipIter = 3;145 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);146 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);147 if (psf->ApTrend == NULL) {148 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");149 return false;150 }151 break;134 // first clip out objects which are too far from the median 135 stats->clipIter = 2; 136 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); 137 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 138 if (psf->ApTrend == NULL) { 139 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); 140 return false; 141 } 142 // apply the fit 143 stats->clipIter = 3; 144 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); 145 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 146 if (psf->ApTrend == NULL) { 147 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); 148 return false; 149 } 150 break; 152 151 case PM_PSF_APTREND_SKYSAT: 153 // first clip out objects which are too far from the median 154 stats->clipIter = 2;155 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);156 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);157 if (psf->ApTrend == NULL) {158 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");159 return false;160 }161 // apply the fit162 stats->clipIter = 2;163 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);164 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);165 if (psf->ApTrend == NULL) {166 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");167 return false;168 }169 // apply the fit170 stats->clipIter = 3;171 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);172 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);173 if (psf->ApTrend == NULL) {174 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");175 return false;176 }177 break;152 // first clip out objects which are too far from the median 153 stats->clipIter = 2; 154 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); 155 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 156 if (psf->ApTrend == NULL) { 157 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); 158 return false; 159 } 160 // apply the fit 161 stats->clipIter = 2; 162 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); 163 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 164 if (psf->ApTrend == NULL) { 165 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); 166 return false; 167 } 168 // apply the fit 169 stats->clipIter = 3; 170 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT); 171 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 172 if (psf->ApTrend == NULL) { 173 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat"); 174 return false; 175 } 176 break; 178 177 case PM_PSF_APTREND_XY_LIN: 179 // first clip out objects which are too far from the median 180 stats->clipIter = 2;181 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);182 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);183 if (psf->ApTrend == NULL) {184 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");185 return false;186 }187 // apply the fit188 stats->clipIter = 3;189 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);190 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);191 if (psf->ApTrend == NULL) {192 psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");193 return false;194 }195 break;178 // first clip out objects which are too far from the median 179 stats->clipIter = 2; 180 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); 181 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 182 if (psf->ApTrend == NULL) { 183 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); 184 return false; 185 } 186 // apply the fit 187 stats->clipIter = 3; 188 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN); 189 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 190 if (psf->ApTrend == NULL) { 191 psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN"); 192 return false; 193 } 194 break; 196 195 case PM_PSF_APTREND_XY_QUAD: 197 // first clip out objects which are too far from the median 198 stats->clipIter = 2;199 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);200 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);201 if (psf->ApTrend == NULL) {202 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");203 return false;204 }205 // apply the fit206 stats->clipIter = 3;207 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);208 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);209 if (psf->ApTrend == NULL) {210 psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");211 return false;212 }213 break;196 // first clip out objects which are too far from the median 197 stats->clipIter = 2; 198 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); 199 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 200 if (psf->ApTrend == NULL) { 201 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); 202 return false; 203 } 204 // apply the fit 205 stats->clipIter = 3; 206 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD); 207 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 208 if (psf->ApTrend == NULL) { 209 psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD"); 210 return false; 211 } 212 break; 214 213 case PM_PSF_APTREND_SKY_XY_LIN: 215 // first clip out objects which are too far from the median 216 stats->clipIter = 2;217 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);218 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);219 if (psf->ApTrend == NULL) {220 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");221 return false;222 }223 // apply the fit224 stats->clipIter = 3;225 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);226 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);227 if (psf->ApTrend == NULL) {228 psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");229 return false;230 }231 break;214 // first clip out objects which are too far from the median 215 stats->clipIter = 2; 216 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); 217 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 218 if (psf->ApTrend == NULL) { 219 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); 220 return false; 221 } 222 // apply the fit 223 stats->clipIter = 3; 224 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN); 225 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 226 if (psf->ApTrend == NULL) { 227 psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin"); 228 return false; 229 } 230 break; 232 231 case PM_PSF_APTREND_SKYSAT_XY_LIN: 233 // first clip out objects which are too far from the median 234 stats->clipIter = 2;235 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);236 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);237 if (psf->ApTrend == NULL) {238 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");239 return false;240 }241 // apply the fit242 stats->clipIter = 3;243 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);244 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);245 if (psf->ApTrend == NULL) {246 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");247 return false;248 }249 // apply the fit250 stats->clipIter = 3;251 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);252 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);253 if (psf->ApTrend == NULL) {254 psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");255 return false;256 }257 break;232 // first clip out objects which are too far from the median 233 stats->clipIter = 2; 234 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); 235 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 236 if (psf->ApTrend == NULL) { 237 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); 238 return false; 239 } 240 // apply the fit 241 stats->clipIter = 3; 242 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); 243 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 244 if (psf->ApTrend == NULL) { 245 psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); 246 return false; 247 } 248 // apply the fit 249 stats->clipIter = 3; 250 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN); 251 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 252 if (psf->ApTrend == NULL) { 253 psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin"); 254 return false; 255 } 256 break; 258 257 case PM_PSF_APTREND_ALL: 259 // first clip out objects which are too far from the median 260 stats->clipIter = 2;261 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);262 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);263 if (psf->ApTrend == NULL) {264 psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");265 return false;266 }267 // fit just SkyBias and clip out objects which are too far from the median 268 stats->clipIter = 2;269 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);270 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);271 if (psf->ApTrend == NULL) {272 psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");273 return false;274 }275 // finally, fit x, y, SkyBias and clip out objects which are too far from the median 276 stats->clipIter = 3;277 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);278 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);279 if (psf->ApTrend == NULL) {280 psError(PSPHOT_ERR_PHOTOM, false, "fitting all");281 return false;282 }283 break;258 // first clip out objects which are too far from the median 259 stats->clipIter = 2; 260 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); 261 psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 262 if (psf->ApTrend == NULL) { 263 psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend"); 264 return false; 265 } 266 // fit just SkyBias and clip out objects which are too far from the median 267 stats->clipIter = 2; 268 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); 269 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 270 if (psf->ApTrend == NULL) { 271 psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias"); 272 return false; 273 } 274 // finally, fit x, y, SkyBias and clip out objects which are too far from the median 275 stats->clipIter = 3; 276 pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL); 277 psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); 278 if (psf->ApTrend == NULL) { 279 psError(PSPHOT_ERR_PHOTOM, false, "fitting all"); 280 return false; 281 } 282 break; 284 283 default: 285 psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName);286 return false;284 psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName); 285 return false; 287 286 } 288 287 … … 295 294 psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV); 296 295 for (int i = 0; i < dMag->n; i++) { 297 if (dMag->data.F64[i] > 0.01) {298 mask->data.U8[i] |= 0x02;299 } 300 if (! mask->data.U8[i]) Nkeep ++;296 if (dMag->data.F64[i] > 0.01) { 297 mask->data.U8[i] |= 0x02; 298 } 299 if (! mask->data.U8[i]) Nkeep ++; 301 300 } 302 301 residStats = psVectorStats (residStats, resid, NULL, mask, 0x03); … … 320 319 321 320 psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot")); 322 psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias, %f skysat (%d of %d used)\n", 323 psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat, Nkeep, Npsf);324 psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n", 325 1e3*psf->ApTrend->coeff[1][0][0][0],326 1e6*psf->ApTrend->coeff[2][0][0][0],327 1e6*psf->ApTrend->coeff[1][1][0][0],328 1e3*psf->ApTrend->coeff[0][1][0][0],329 1e6*psf->ApTrend->coeff[0][2][0][0]);321 psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias, %f skysat (%d of %d used)\n", 322 psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat, Nkeep, Npsf); 323 psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n", 324 1e3*psf->ApTrend->coeff[1][0][0][0], 325 1e6*psf->ApTrend->coeff[2][0][0][0], 326 1e6*psf->ApTrend->coeff[1][1][0][0], 327 1e3*psf->ApTrend->coeff[0][1][0][0], 328 1e6*psf->ApTrend->coeff[0][2][0][0]); 330 329 331 330 psFree (mask);
Note:
See TracChangeset
for help on using the changeset viewer.
