Changeset 9595 for trunk/psphot/src/psphotEnsemblePSF.c
- Timestamp:
- Oct 16, 2006, 4:24:03 PM (20 years ago)
- File:
-
- 1 edited
-
trunk/psphot/src/psphotEnsemblePSF.c (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot/src/psphotEnsemblePSF.c
r9576 r9595 4 4 // 2006.02.07 : no leaks! 5 5 // fit all reasonable sources with the linear PSF model 6 bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { 6 bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { 7 7 8 8 bool status; … … 39 39 40 40 const bool CONSTANT_PHOTOMETRIC_WEIGHTS = 41 41 psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); 42 42 if (!status) { 43 43 psAbort(PS_FILE_LINE, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); 44 44 } 45 45 46 46 for (int i = 0; i < sources->n; i++) { 47 48 49 50 51 if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue; 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); 75 76 77 78 79 80 81 // ix,iy must land on inSource->pixels 82 83 84 85 86 87 88 psTrace ("psphot", 5, "peak coord: %f %f -> %f %f\n", 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 // build the model image 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 47 pmSource *inSource = sources->data[i]; 48 49 // skip non-astronomical objects (very likely defects) 50 // XXX EAM : should we try these anyway? 51 if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue; 52 if (inSource->type == PM_SOURCE_TYPE_SATURATED) continue; 53 if (final) { 54 if (inSource->mode & PM_SOURCE_MODE_SUBTRACTED) continue; 55 } else { 56 if (inSource->mode & PM_SOURCE_MODE_BLEND) continue; 57 } 58 59 if (inSource->moments->x < AnalysisRegion.x0) continue; 60 if (inSource->moments->y < AnalysisRegion.y0) continue; 61 if (inSource->moments->x > AnalysisRegion.x1) continue; 62 if (inSource->moments->y > AnalysisRegion.y1) continue; 63 64 pmSource *otSource = pmSourceAlloc (); 65 66 // really saturated stars should be re-measured for a better centroid 67 // XXX EAM : move this to a 'clear satstar function' 68 if (inSource->mode & PM_SOURCE_MODE_SATSTAR) { 69 status = pmSourceMoments (inSource, INNER_RADIUS); 70 } 71 72 // XXX EAM : add option to use EXT or PSF form 73 // use the source moments, etc to guess basic model parameters 74 pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); 75 if (inSource->mode & PM_SOURCE_MODE_SATSTAR) { 76 modelEXT->params->data.F32[PM_PAR_XPOS] = inSource->moments->x; 77 modelEXT->params->data.F32[PM_PAR_YPOS] = inSource->moments->y; 78 } else { 79 // peak-up on peak (for non-sat objects) 80 81 // ix,iy must land on inSource->pixels 82 int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x)); 83 int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y)); 84 85 psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy); 86 psPlane min = psImageBicubeMin (bicube); 87 88 psTrace ("psphot", 5, "peak coord: %f %f -> %f %f\n", 89 modelEXT->params->data.F32[PM_PAR_XPOS], modelEXT->params->data.F32[PM_PAR_YPOS], min.x + ix, min.y + iy); 90 91 // if min point is too deviant, keep the old value 92 if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) { 93 modelEXT->params->data.F32[PM_PAR_XPOS] = min.x + ix; 94 modelEXT->params->data.F32[PM_PAR_YPOS] = min.y + iy; 95 } 96 psFree (bicube); 97 } 98 99 // set PSF parameters for this model 100 pmModel *model = pmModelFromPSF (modelEXT, psf); 101 psFree (modelEXT); 102 103 // save the original coords 104 x = model->params->data.F32[PM_PAR_XPOS]; 105 y = model->params->data.F32[PM_PAR_YPOS]; 106 107 // set the fit radius based on the object flux limit and the model 108 psphotCheckRadiusPSF (readout, inSource, model); 109 110 // make temporary copies of the image pixels and mask 111 otSource->mask = psImageCopy (NULL, inSource->mask, PS_TYPE_U8); 112 otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32); 113 otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32); 114 115 // build the model image 116 psImage *flux = otSource->pixels; 117 psImage *mask = otSource->mask; 118 119 // set model to unit peak, zero sky (we assume sky is subtracted) 120 model->params->data.F32[PM_PAR_SKY] = 0.0; 121 model->params->data.F32[PM_PAR_I0] = 1.0; 122 123 // fill in the model pixel values 124 psImageInit (flux, 0.0); 125 psImageKeepCircle (mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK); 126 pmModelAdd (flux, mask, model, false, false); 127 128 // calculate nDOF (nPix - 1) 129 // int Nmaskpix = psImageCountPixelMask (mask, allArray, PM_MASK_SAT); 130 // model->nDOF = mask->numCols*mask->numRows - Nmaskpix - 1; 131 132 // save source in list 133 otSource->modelPSF = model; 134 index->data.U32[models->n] = i; 135 psArrayAdd (models, 100, otSource); 136 psFree (otSource); 137 137 } 138 138 psLogMsg ("psphot.emsemble", 4, "built models: %f (%ld objects)\n", psTimerMark ("psphot"), sources->n); 139 139 140 140 // fill out the sparse matrix 141 141 psSparse *sparse = psSparseAlloc (models->n, 100); … … 144 144 145 145 for (int i = 0; i < models->n; i++) { 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 psSparseMatrixElement (sparse, j, i, f / r); 173 146 int N = index->data.U32[i]; 147 pmSource *Fi = sources->data[N]; 148 pmSource *Mi = models->data[i]; 149 150 // scale by diagonal element (auto-cross-product) 151 r = pmSourceCrossProduct (Mi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS); 152 weight->data.F32[i] = r; 153 154 psSparseMatrixElement (sparse, i, i, 1.0); 155 156 // find the image x model value 157 f = pmSourceCrossProduct (Fi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS); 158 psSparseVectorElement (sparse, i, f / r); 159 160 // loop over all other stars following this one 161 for (int j = i + 1; j < models->n; j++) { 162 pmSource *Mj = models->data[j]; 163 164 // skip over disjoint source images, break after last possible overlap 165 if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break; 166 if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue; 167 if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue; 168 if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue; 169 170 // got an overlap; calculate cross-product and add to output array 171 f = pmSourceCrossProduct (Mi, Mj, CONSTANT_PHOTOMETRIC_WEIGHTS); 172 psSparseMatrixElement (sparse, j, i, f / r); 173 } 174 174 } 175 175 psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem); … … 186 186 // adjust models, set sources and subtract 187 187 for (int i = 0; i < models->n; i++) { 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 188 int N = index->data.U32[i]; 189 pmSource *Fi = sources->data[N]; 190 pmSource *Mi = models->data[i]; 191 192 // if we already have a PSF model, free it. 193 psFree (Fi->modelPSF); 194 195 // need to increment counter so we can free models here and sources above 196 Fi->modelPSF = psMemIncrRefCounter (Mi->modelPSF); 197 198 // assign linearly-fitted normalization 199 if (isnan(norm->data.F32[i])) { 200 psAbort ("psphot", "ensemble source is nan"); 201 } 202 Fi->modelPSF->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; 203 Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i])); 204 // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it 205 206 // subtract object 207 pmModelSub (Fi->pixels, Fi->mask, Fi->modelPSF, false, false); 208 Fi->mode |= PM_SOURCE_MODE_SUBTRACTED; 209 if (!final) Fi->mode |= PM_SOURCE_MODE_TEMPSUB; 210 210 } 211 211 212 212 // measure chisq for each source 213 213 for (int i = 0; final && (i < models->n); i++) { 214 215 216 217 218 219 220 221 222 223 psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "AND",NOT_U8(PM_MASK_MARK));214 int N = index->data.U32[i]; 215 pmSource *Fi = sources->data[N]; 216 pmModel *model = Fi->modelPSF; 217 218 x = model->params->data.F32[PM_PAR_XPOS]; 219 y = model->params->data.F32[PM_PAR_YPOS]; 220 221 psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK); 222 pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight); 223 psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK)); 224 224 } 225 225 … … 239 239 int Npix = 0; 240 240 for (int j = 0; j < image->numRows; j++) { 241 242 243 244 245 246 241 for (int i = 0; i < image->numCols; i++) { 242 if (mask->data.U8[j][i]) continue; 243 if (weight->data.F32[j][i] <= 0) continue; 244 dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i]; 245 Npix ++; 246 } 247 247 } 248 248 model->nDOF = Npix - 1;
Note:
See TracChangeset
for help on using the changeset viewer.
