Changeset 8046
- Timestamp:
- Aug 1, 2006, 4:15:32 PM (20 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
-
psastro/src/psastroAstromGuess.c (modified) (6 diffs)
-
psphot/src/psphotChoosePSF.c (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psastro/src/psastroAstromGuess.c
r7332 r8046 4 4 // XXX WCS keywords in the headers corresponding to the chips. 5 5 // XXX this function is both converting the header WCS astrometry terms to the fpa terms 6 // and applying the astrometry to the detected objects. 6 // and applying the astrometry to the detected objects. 7 7 bool psastroAstromGuess (pmConfig *config) { 8 8 … … 10 10 bool status = false; 11 11 bool isMosaic = false; 12 double RAmin , RAmax, RAminSky, RAmaxSky;13 double DECmin , DECmax;12 double RAmin = NAN, RAmax = NAN, RAminSky = NAN, RAmaxSky = NAN; 13 double DECmin = NAN, DECmax = NAN; 14 14 15 15 pmChip *chip = NULL; … … 23 23 psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO"); 24 24 if (!recipe) { 25 psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n");26 exit(EXIT_FAILURE);25 psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n"); 26 exit(EXIT_FAILURE); 27 27 } 28 28 … … 30 30 pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); 31 31 if (!input) { 32 psErrorStackPrint(stderr, "Can't find input data!\n");33 exit(EXIT_FAILURE);32 psErrorStackPrint(stderr, "Can't find input data!\n"); 33 exit(EXIT_FAILURE); 34 34 } 35 35 … … 44 44 psTrace (__func__, 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); 45 45 if (!chip->process || !chip->file_exists) { continue; } 46 46 47 47 // read WCS data from the corresponding header 48 48 pmHDU *hdu = pmFPAviewThisHDU (view, fpa); 49 49 50 pmAstromReadWCS (fpa, chip, hdu->header, plateScale, isMosaic);51 if (newFPA) {52 newFPA = false;53 RAminSky = fpa->projection->R - M_PI;54 RAmaxSky = fpa->projection->R + M_PI;55 RAmin = RAmax = fpa->projection->R;56 DECmin = DECmax = fpa->projection->D;57 }50 pmAstromReadWCS (fpa, chip, hdu->header, plateScale, isMosaic); 51 if (newFPA) { 52 newFPA = false; 53 RAminSky = fpa->projection->R - M_PI; 54 RAmaxSky = fpa->projection->R + M_PI; 55 RAmin = RAmax = fpa->projection->R; 56 DECmin = DECmax = fpa->projection->D; 57 } 58 58 59 // apply the new WCS guess data to all of the data in the readouts60 // XXX should this go into a different function? this would separate WCS interpretation from application61 while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {59 // apply the new WCS guess data to all of the data in the readouts 60 // XXX should this go into a different function? this would separate WCS interpretation from application 61 while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { 62 62 psTrace (__func__, 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); 63 63 if (!cell->process || !cell->file_exists) { continue; } 64 64 65 // process each of the readouts66 while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {67 if (! readout->data_exists) { continue; }65 // process each of the readouts 66 while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { 67 if (! readout->data_exists) { continue; } 68 68 69 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");70 if (rawstars == NULL) { continue; }69 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); 70 if (rawstars == NULL) { continue; } 71 71 72 for (int i = 0; i < rawstars->n; i++) { 73 pmAstromObj *raw = rawstars->data[i]; 74 75 psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); 76 psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0); 77 p_psDeproject (raw->sky, raw->TP, fpa->projection); 72 for (int i = 0; i < rawstars->n; i++) { 73 pmAstromObj *raw = rawstars->data[i]; 78 74 79 if (i < 0) { 80 fprintf (stderr, "up: %f,%f -> %f,%f -> %f,%f -> %f,%f\n", 81 raw->chip->x, raw->chip->y, 82 raw->FP->x, raw->FP->y, 83 raw->TP->x, raw->TP->y, 84 raw->sky->r, raw->sky->d); 75 psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); 76 psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0); 77 p_psDeproject (raw->sky, raw->TP, fpa->projection); 85 78 86 psPlane *fp = psPlaneAlloc(); 87 psPlane *tp = psPlaneAlloc(); 88 psPlane *ch = psPlaneAlloc(); 79 if (i < 0) { 80 fprintf (stderr, "up: %f,%f -> %f,%f -> %f,%f -> %f,%f\n", 81 raw->chip->x, raw->chip->y, 82 raw->FP->x, raw->FP->y, 83 raw->TP->x, raw->TP->y, 84 raw->sky->r, raw->sky->d); 89 85 90 p_psProject (tp, raw->sky, fpa->projection);91 psPlaneDistortApply (fp, fpa->fromTangentPlane, tp, 0.0, 0.0);92 psPlaneTransformApply (ch, chip->fromFPA, fp);86 psPlane *fp = psPlaneAlloc(); 87 psPlane *tp = psPlaneAlloc(); 88 psPlane *ch = psPlaneAlloc(); 93 89 94 fprintf (stderr, "dn: %f,%f <- %f,%f <- %f,%f <- %f,%f\n", 95 ch->x, ch->y, 96 fp->x, fp->y, 97 tp->x, tp->y, 98 raw->sky->r, raw->sky->d); 99 } 90 p_psProject (tp, raw->sky, fpa->projection); 91 psPlaneDistortApply (fp, fpa->fromTangentPlane, tp, 0.0, 0.0); 92 psPlaneTransformApply (ch, chip->fromFPA, fp); 100 93 101 // rationalize ra to sky range centered on boresite 102 while (raw->sky->r < RAminSky) raw->sky->r += M_PI; 103 while (raw->sky->r > RAmaxSky) raw->sky->r -= M_PI; 94 fprintf (stderr, "dn: %f,%f <- %f,%f <- %f,%f <- %f,%f\n", 95 ch->x, ch->y, 96 fp->x, fp->y, 97 tp->x, tp->y, 98 raw->sky->r, raw->sky->d); 99 } 104 100 105 RAmin = PS_MIN (raw->sky->r, RAmin); 106 RAmax = PS_MAX (raw->sky->r, RAmax); 101 // rationalize ra to sky range centered on boresite 102 while (raw->sky->r < RAminSky) raw->sky->r += M_PI; 103 while (raw->sky->r > RAmaxSky) raw->sky->r -= M_PI; 107 104 108 DECmin = PS_MIN (raw->sky->d, DECmin); 109 DECmax = PS_MAX (raw->sky->d, DECmax); 110 } 111 } 112 } 105 RAmin = PS_MIN (raw->sky->r, RAmin); 106 RAmax = PS_MAX (raw->sky->r, RAmax); 107 108 DECmin = PS_MIN (raw->sky->d, DECmin); 109 DECmax = PS_MAX (raw->sky->d, DECmax); 110 } 111 } 112 } 113 113 } 114 114 … … 119 119 psMetadataAdd (recipe, PS_LIST_TAIL, "DEC_MIN", PS_DATA_F32 | PS_META_REPLACE, "", DECmin); 120 120 psMetadataAdd (recipe, PS_LIST_TAIL, "DEC_MAX", PS_DATA_F32 | PS_META_REPLACE, "", DECmax); 121 121 122 122 psFree (view); 123 123 return true; -
trunk/psphot/src/psphotChoosePSF.c
r7758 r8046 2 2 3 3 // try PSF models and select best option 4 pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) { 4 pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) { 5 5 6 bool status;7 char *modelName;8 pmPSF *psf = NULL;9 pmPSFtry *try = NULL;10 psArray *stars = NULL;6 bool status; 7 char *modelName; 8 pmPSF *psf = NULL; 9 pmPSFtry *try = NULL; 10 psArray *stars = NULL; 11 11 12 12 psTimerStart ("psphot"); … … 32 32 // select the candidate PSF stars (pointers to original sources) 33 33 for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) { 34 pmSource *source = sources->data[i];35 if (source->mode & PM_SOURCE_MODE_PSFSTAR) psArrayAdd (stars, 200, source);34 pmSource *source = sources->data[i]; 35 if (source->mode & PM_SOURCE_MODE_PSFSTAR) psArrayAdd (stars, 200, source); 36 36 } 37 37 psLogMsg ("psphot.pspsf", 4, "selected candidate %d PSF objects\n", stars->n); … … 47 47 48 48 if (mdi->type == PS_DATA_STRING) { 49 list = psListAlloc(NULL);50 psListAdd (list, PS_LIST_HEAD, mdi);49 list = psListAlloc(NULL); 50 psListAdd (list, PS_LIST_HEAD, mdi); 51 51 } else { 52 if (mdi->type != PS_DATA_METADATA_MULTI) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");53 list = psMemIncrRefCounter(mdi->data.list);52 if (mdi->type != PS_DATA_METADATA_MULTI) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection"); 53 list = psMemIncrRefCounter(mdi->data.list); 54 54 } 55 55 … … 59 59 60 60 // try each model option listed in config 61 psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); 62 for (int i = 0; i < models->n; i++) { 63 psMetadataItem *item = psListGetAndIncrement (iter);64 modelName = item->data.V;65 models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS);61 psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); 62 for (int i = 0; i < models->n; i++) { 63 psMetadataItem *item = psListGetAndIncrement (iter); 64 modelName = item->data.V; 65 models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS); 66 66 } 67 67 psFree (iter); … … 75 75 float bestM = try->psf->dApResid; 76 76 for (int i = 1; i < models->n; i++) { 77 try = models->data[i];78 float M = try->psf->dApResid;79 if (M < bestM) {80 bestM = M;81 bestN = i;82 }77 try = models->data[i]; 78 float M = try->psf->dApResid; 79 if (M < bestM) { 80 bestM = M; 81 bestN = i; 82 } 83 83 } 84 84 85 85 // use the best model: 86 86 try = models->data[bestN]; … … 88 88 // unset the PSFSTAR flag for stars not used for PSF model 89 89 for (int i = 0; i < try->sources->n; i++) { 90 pmSource *source = try->sources->data[i];91 if (try->mask->data.U8[i] & PSFTRY_MASK_EXT_FAIL) {92 source->mode &= ~PM_SOURCE_MODE_PSFSTAR;93 }90 pmSource *source = try->sources->data[i]; 91 if (try->mask->data.U8[i] & PSFTRY_MASK_EXT_FAIL) { 92 source->mode &= ~PM_SOURCE_MODE_PSFSTAR; 93 } 94 94 } 95 95 … … 130 130 // get the model full-width at half-max 131 131 psF64 FWHM_X = 2*modelRadiusFunc (modelPSF->params, 0.5); 132 132 133 133 shape.sx = modelPSF->params->data.F32[4]; 134 134 shape.sy = modelPSF->params->data.F32[5]; … … 141 141 psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y); 142 142 psMetadataAdd (recipe, PS_LIST_TAIL, "ANGLE", PS_DATA_F32 | PS_META_REPLACE, "PSF angle", axes.theta); 143 143 144 144 psFree (modelEXT); 145 145 psFree (modelPSF); … … 150 150 bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources) { 151 151 152 // without the PSF model, we can only rely on the source->moments 152 // without the PSF model, we can only rely on the source->moments 153 153 // as a measure of the FWHM values 154 154 155 double FWHM_X , FWHM_Y, FWHM_T;156 int FWHM_N ;155 double FWHM_X = 0.0, FWHM_Y = 0.0, FWHM_T = 0.0; 156 int FWHM_N = 0; 157 157 158 158 psEllipseMoments moments; 159 159 psEllipseAxes axes; 160 160 161 FWHM_N = 0; 162 FWHM_X = FWHM_Y = 0.0; 161 for (int i = 0; i < sources->n; i++) { 162 pmSource *source = sources->data[i]; 163 if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; 163 164 164 for (int i = 0; i < sources->n; i++) { 165 pmSource *source = sources->data[i]; 166 if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; 167 168 // moments->Sx,Sy are roughly sigma_x,y 169 moments.x2 = PS_SQR(source->moments->Sx); 170 moments.y2 = PS_SQR(source->moments->Sy); 171 moments.xy = source->moments->Sxy; 165 // moments->Sx,Sy are roughly sigma_x,y 166 moments.x2 = PS_SQR(source->moments->Sx); 167 moments.y2 = PS_SQR(source->moments->Sy); 168 moments.xy = source->moments->Sxy; 172 169 173 axes = psEllipseMomentsToAxes (moments);170 axes = psEllipseMomentsToAxes (moments); 174 171 175 FWHM_X += axes.major * 2.35;176 FWHM_Y += axes.minor * 2.35;177 FWHM_T += axes.theta;178 FWHM_N ++;172 FWHM_X += axes.major * 2.35; 173 FWHM_Y += axes.minor * 2.35; 174 FWHM_T += axes.theta; 175 FWHM_N ++; 179 176 } 180 177 … … 186 183 psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y); 187 184 psMetadataAdd (recipe, PS_LIST_TAIL, "ANGLE", PS_DATA_F32 | PS_META_REPLACE, "PSF angle", FWHM_T); 188 185 189 186 return true; 190 187 }
Note:
See TracChangeset
for help on using the changeset viewer.
