Changeset 28440 for branches/eam_branches/ipp-20100621/psphot
- Timestamp:
- Jun 23, 2010, 2:36:55 PM (16 years ago)
- Location:
- branches/eam_branches/ipp-20100621/psphot/src
- Files:
-
- 6 edited
-
psphotFitSourcesLinear.c (modified) (4 diffs)
-
psphotFitSourcesLinearStack.c (modified) (5 diffs)
-
psphotImageLoop.c (modified) (3 diffs)
-
psphotRadiusChecks.c (modified) (6 diffs)
-
psphotRoughClass.c (modified) (5 diffs)
-
psphotSourceStats.c (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinear.c
r28399 r28440 171 171 172 172 // diagonal elements of the sparse matrix (auto-cross-product) 173 f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor );173 f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); 174 174 psSparseMatrixElement (sparse, i, i, f); 175 175 176 176 // the formal error depends on the weighting scheme 177 177 if (CONSTANT_PHOTOMETRIC_WEIGHTS) { 178 float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor );178 float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal); 179 179 errors->data.F32[i] = 1.0 / sqrt(var); 180 180 } else { … … 184 184 185 185 // find the image x model value 186 f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor );186 f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); 187 187 psSparseVectorElement (sparse, i, f); 188 188 … … 190 190 switch (SKY_FIT_ORDER) { 191 191 case 1: 192 f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor );192 f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); 193 193 psSparseBorderElementB (border, i, 1, f); 194 f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor );194 f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); 195 195 psSparseBorderElementB (border, i, 2, f); 196 196 197 197 case 0: 198 f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor );198 f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); 199 199 psSparseBorderElementB (border, i, 0, f); 200 200 break; … … 216 216 217 217 // got an overlap; calculate cross-product and add to output array 218 f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor );218 f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); 219 219 psSparseMatrixElement (sparse, j, i, f); 220 220 } -
branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinearStack.c
r28013 r28440 43 43 for (int i = 0; i < objects->n; i++) { 44 44 pmPhotObj *object = objects->data[i]; 45 if (!object) continue;46 if (!object->sources) continue;45 if (!object) continue; 46 if (!object->sources) continue; 47 47 48 // XXX check an element of the group to see if we should use it49 // if (!object->flags & PM_PHOT_OBJ_BAD) continue;48 // XXX check an element of the group to see if we should use it 49 // if (!object->flags & PM_PHOT_OBJ_BAD) continue; 50 50 51 for (int j = 0; j < object->sources->n; j++) {52 pmSource *source = object->sources->data[j];53 if (!source) continue;51 for (int j = 0; j < object->sources->n; j++) { 52 pmSource *source = object->sources->data[j]; 53 if (!source) continue; 54 54 55 // turn this bit off and turn it on again if we keep this source56 source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;55 // turn this bit off and turn it on again if we keep this source 56 source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT; 57 57 58 // generate model for sources without, or skip if we can't59 if (!source->modelFlux) {58 // generate model for sources without, or skip if we can't 59 if (!source->modelFlux) { 60 60 if (!pmSourceCacheModel (source, maskVal)) continue; 61 }61 } 62 62 63 source->mode |= PM_SOURCE_MODE_LINEAR_FIT;64 psArrayAdd (fitSources, 100, source);65 }63 source->mode |= PM_SOURCE_MODE_LINEAR_FIT; 64 psArrayAdd (fitSources, 100, source); 65 } 66 66 } 67 67 psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n); … … 85 85 86 86 // diagonal elements of the sparse matrix (auto-cross-product) 87 f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR );87 f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal); 88 88 psSparseMatrixElement (sparse, i, i, f); 89 89 90 90 // the formal error depends on the weighting scheme 91 91 if (CONSTANT_PHOTOMETRIC_WEIGHTS) { 92 float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR );92 float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR, maskVal); 93 93 errors->data.F32[i] = 1.0 / sqrt(var); 94 94 } else { … … 97 97 98 98 // find the image x model value 99 f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR );99 f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal); 100 100 psSparseVectorElement (sparse, i, f); 101 101 … … 104 104 pmSource *SRCj = fitSources->data[j]; 105 105 106 // we only need to generate dot terms for source on the same image107 if (SRCj->imageID != SRCi->imageID) { continue; }106 // we only need to generate dot terms for source on the same image 107 if (SRCj->imageID != SRCi->imageID) { continue; } 108 108 109 109 // skip over disjoint source images, break after last possible overlap … … 114 114 115 115 // got an overlap; calculate cross-product and add to output array 116 f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR );116 f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal); 117 117 psSparseMatrixElement (sparse, j, i, f); 118 118 } -
branches/eam_branches/ipp-20100621/psphot/src/psphotImageLoop.c
r27657 r28440 46 46 if (!psphotMosaicChip(config, view, "PSPHOT.INPUT", "PSPHOT.LOAD")) ESCAPE ("Unable to mosaic chip."); 47 47 48 // Read WCS if easy. 49 // XXX Since we're mosaicking cells, we ignore the case where the WCS is defined for a cell. 50 { 51 pmChip *inChip = pmFPAviewThisChip(view, input->fpa); // Mosaicked chip 52 pmHDU *hduLow = pmHDUGetLowest(input->fpa, inChip, NULL); 53 if (hduLow && !pmAstromReadWCS(input->fpa, inChip, hduLow->header, 1.0)) { 54 psWarning("Unable to read WCS astrometry from header."); 55 psErrorClear(); 56 pmHDU *hduHigh = pmHDUGetHighest(input->fpa, inChip, NULL); 57 if (hduHigh && hduHigh != hduLow && 58 !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) { 59 psWarning("Unable to read WCS astrometry from primary header."); 60 psErrorClear(); 61 } 62 } 63 } 64 48 65 // try to load other supporting data (PSF, SRC, etc). 49 66 // do not re-load the following three files … … 67 84 68 85 // Update the header 69 pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); 70 if (hdu && hdu != lastHDU) { 71 psphotVersionHeaderFull(hdu->header); 72 lastHDU = hdu; 86 { 87 pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); 88 if (hdu && hdu != lastHDU) { 89 psphotVersionHeaderFull(hdu->header); 90 lastHDU = hdu; 91 } 73 92 } 74 93 75 // if an external mask is supplied, ensure that NAN pixels are also masked76 if (readout->mask) {77 psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels78 if (!pmReadoutMaskNonfinite(readout, maskSat)) {79 psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");80 psFree(view);81 return false;82 }83 }94 // if an external mask is supplied, ensure that NAN pixels are also masked 95 if (readout->mask) { 96 psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels 97 if (!pmReadoutMaskNonfinite(readout, maskSat)) { 98 psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); 99 psFree(view); 100 return false; 101 } 102 } 84 103 85 104 // run the actual photometry analysis on this chip/cell/readout … … 91 110 } 92 111 93 // drop all versions of the internal files94 status = true;95 status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");96 status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");97 status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");98 if (!status) {99 psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files");100 psFree (view);101 return false;102 }112 // drop all versions of the internal files 113 status = true; 114 status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); 115 status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); 116 status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); 117 if (!status) { 118 psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); 119 psFree (view); 120 return false; 121 } 103 122 } 104 123 // save output which is saved at the chip level -
branches/eam_branches/ipp-20100621/psphot/src/psphotRadiusChecks.c
r26894 r28440 4 4 static float PSF_FIT_NSIGMA; 5 5 static float PSF_FIT_PADDING; 6 static float PSF_APERTURE = 0; // radius to use in PSF aperture mags7 static float PSF_FIT_RADIUS = 0; // radius to use in fitting (ignored if <= 0,8 // and a per-object radius is calculated)6 static float PSF_APERTURE = 0; // radius to use in PSF aperture mags 7 static float PSF_FIT_RADIUS = 0; // radius to use in fitting (ignored if <= 0, 8 // and a per-object radius is calculated) 9 9 10 10 bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) { … … 17 17 PSF_FIT_RADIUS = psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS"); 18 18 if (!status) { 19 PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");19 PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); 20 20 } 21 21 22 22 PSF_APERTURE = psMetadataLookupF32(&status, analysis, "PSF_APERTURE"); 23 23 if (!status) { 24 PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); 24 PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); 25 } 26 27 // The PSF_FIT_RADIUS and PSF_APERTURE may not be set if the PSF was loaded and not chosen 28 29 if (PSF_FIT_RADIUS == 0.0) { 30 float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA"); 31 if (!status) { 32 gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); 33 } 34 float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE"); 35 PSF_FIT_RADIUS = (int)(fitScale*gaussSigma); 36 } 37 38 if (PSF_APERTURE == 0.0) { 39 float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA"); 40 if (!status) { 41 gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); 42 } 43 float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE"); 44 PSF_APERTURE = (int)(apScale*gaussSigma); 25 45 } 26 46 … … 38 58 // set the fit radius based on the object flux limit and the model 39 59 float radiusFit = PSF_FIT_RADIUS; 40 if (radiusFit <= 0) { // use fixed radius41 if (moments == NULL) {42 radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);43 } else {44 radiusFit = model->modelRadius(model->params, 1.0);45 }46 model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);60 if (radiusFit <= 0) { // use fixed radius 61 if (moments == NULL) { 62 radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); 63 } else { 64 radiusFit = model->modelRadius(model->params, 1.0); 65 } 66 model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); 47 67 } else { 48 model->fitRadius = radiusFit;68 model->fitRadius = radiusFit; 49 69 } 50 70 if (isnan(model->fitRadius)) psAbort("error in radius"); 51 71 52 72 if (source->mode & PM_SOURCE_MODE_SATSTAR) { 53 model->fitRadius *= 2;73 model->fitRadius *= 2; 54 74 } 55 75 … … 73 93 // set the fit radius based on the object flux limit and the model 74 94 float radiusFit = PSF_FIT_RADIUS; 75 if (radiusFit <= 0) { // use fixed radius76 if (moments == NULL) {77 radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);78 } else {79 radiusFit = model->modelRadius(model->params, 1.0);80 }81 model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);95 if (radiusFit <= 0) { // use fixed radius 96 if (moments == NULL) { 97 radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); 98 } else { 99 radiusFit = model->modelRadius(model->params, 1.0); 100 } 101 model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); 82 102 } else { 83 model->fitRadius = radiusFit;103 model->fitRadius = radiusFit; 84 104 } 85 105 if (isnan(model->fitRadius)) psAbort("error in radius"); … … 89 109 90 110 if (source->mode & PM_SOURCE_MODE_SATSTAR) { 91 model->fitRadius *= 2;111 model->fitRadius *= 2; 92 112 } 93 113 … … 134 154 float radius = 0.0; 135 155 for (int j = 0; j < footprint->spans->n; j++) { 136 pmSpan *span = footprint->spans->data[j];137 138 float dY = span->y - peak->yf;139 float dX0 = span->x0 - peak->xf;140 float dX1 = span->x1 - peak->xf;141 142 radius = PS_MAX (radius, hypot(dY, dX0));143 radius = PS_MAX (radius, hypot(dY, dX1));156 pmSpan *span = footprint->spans->data[j]; 157 158 float dY = span->y - peak->yf; 159 float dX0 = span->x0 - peak->xf; 160 float dX1 = span->x1 - peak->xf; 161 162 radius = PS_MAX (radius, hypot(dY, dX0)); 163 radius = PS_MAX (radius, hypot(dY, dX1)); 144 164 } 145 165 -
branches/eam_branches/ipp-20100621/psphot/src/psphotRoughClass.c
r28013 r28440 25 25 // loop over the available readouts 26 26 for (int i = 0; i < num; i++) { 27 if (i == chisqNum) continue; // skip chisq image28 if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) {27 if (i == chisqNum) continue; // skip chisq image 28 if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) { 29 29 psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for %s entry %d", filerule, i); 30 return false;31 }30 return false; 31 } 32 32 } 33 33 return true; … … 50 50 bool havePSF = false; 51 51 if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { 52 havePSF = true;52 havePSF = true; 53 53 } 54 54 … … 60 60 61 61 if (!sources->n) { 62 psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification");63 return true;62 psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification"); 63 return true; 64 64 } 65 65 … … 78 78 psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n", 79 79 region->x0, region->y0, region->x1, region->y1); 80 81 // If in doubt, it's a PSF 82 for (int i = 0; i < sources->n; i++) { 83 pmSource *source = sources->data[i]; // Source of interest 84 if (!source || !source->peak) { 85 continue; 86 } 87 if (source->peak->x < region->x0) continue; 88 if (source->peak->x >= region->x1) continue; 89 if (source->peak->y < region->y0) continue; 90 if (source->peak->y >= region->y1) continue; 91 source->type = PM_SOURCE_TYPE_STAR; 92 } 80 93 psFree (region); 81 94 continue; … … 124 137 // XXX why not save the psfClump as a PTR? 125 138 126 float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");127 float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");139 float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); 140 float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX"); 128 141 129 float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE");130 if (!status) {131 PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");132 psAssert (status, "missing PSF_CLUMP_GRID_SCALE");133 }134 float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX");135 if (!status) {136 MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");137 psAssert (status, "missing MOMENTS_SX_MAX");138 }139 float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX");140 if (!status) {141 MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");142 psAssert (status, "missing MOMENTS_SY_MAX");143 }142 float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE"); 143 if (!status) { 144 PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE"); 145 psAssert (status, "missing PSF_CLUMP_GRID_SCALE"); 146 } 147 float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX"); 148 if (!status) { 149 MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX"); 150 psAssert (status, "missing MOMENTS_SX_MAX"); 151 } 152 float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX"); 153 if (!status) { 154 MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX"); 155 psAssert (status, "missing MOMENTS_SY_MAX"); 156 } 144 157 145 158 psfClump = pmSourcePSFClump (NULL, region, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX); -
branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c
r28435 r28440 489 489 psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); 490 490 491 #if 0 492 // Modifying clump parameters without restoring! 491 493 psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); 492 494 psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, "PSF.CLUMP.REGION.000"); … … 500 502 psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); 501 503 psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); 502 503 504 if (pmVisualTestLevel("psphot.moments.full", 2)) { 504 505 psphotVisualPlotMoments (recipe, analysis, sources); 505 506 } 507 #endif 506 508 507 509 Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
Note:
See TracChangeset
for help on using the changeset viewer.
