Changeset 15016 for trunk/psModules/src/objects/pmPSFtry.c
- Timestamp:
- Sep 25, 2007, 12:05:05 PM (19 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/objects/pmPSFtry.c (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/pmPSFtry.c
r15002 r15016 5 5 * @author EAM, IfA 6 6 * 7 * @version $Revision: 1.4 7$ $Name: not supported by cvs2svn $8 * @date $Date: 2007-09-25 04:17:29$7 * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2007-09-25 22:05:05 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 35 35 #include "pmSourcePhotometry.h" 36 36 37 int testSaveImage (psMetadata *header, psImage *image, char *filename) {38 39 psFits *fits = psFitsOpen (filename, "w");40 psFitsWriteImage (fits, NULL, image, 0, NULL);41 psFitsClose (fits);42 return (TRUE);43 }44 45 bool pmPSFParamTrend (pmPSF *psf, int scale, float *errorFloor, float *errorTotal, psVector *mask, psVector *x, psVector *y, psVector *e0, psVector *e1, psVector *e2, psVector *dz);46 bool pmPSFErrorFloor (float *errorFloor, psVector *dMag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup);47 48 37 // ******** pmPSFtry functions ************************************************** 49 38 // * pmPSFtry holds a single pmPSF model test, with the input sources, the freely … … 363 352 psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); 364 353 psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); 365 psVector *flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);366 367 psVector *dz = NULL;368 if (psf->poissonErrorsParams) {369 dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);370 }371 354 372 355 // construct the x,y terms … … 378 361 x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; 379 362 y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; 380 flux->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_I0]; 381 382 // weight by the error on the source flux 383 if (dz) { 384 // dz->data.F32[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0] / source->modelEXT->params->data.F32[PM_PAR_I0]; 385 dz->data.F32[i] = 1.0; 386 } 387 } 388 389 // we are doing a robust fit. after each pass, we drop points which are 390 // more deviant than three sigma. 391 // mask is currently updated for each pass. this is inconsistent? 392 393 // The shape parameters (SXX, SXY, SYY) are strongly coupled. We have to handle them very 394 // carefully. First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for 395 // each source and fit this set of parameters. These values are less tightly coupled, but 396 // are still inter-related. The fitted values do a good job of constraining the major axis 397 // and the position angle, but the minor axis is weakly measured. When we apply the PSF 398 // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape 399 // parameters, with the constraint that the minor axis must be greater than a minimum 400 // threshold. 401 402 // convert the measured source shape paramters to polarization terms 403 psVector *e0 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); 404 psVector *e1 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); 405 psVector *e2 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); 406 for (int i = 0; i < psfTry->sources->n; i++) { 407 pmSource *source = psfTry->sources->data[i]; 408 if (source->modelEXT == NULL) continue; 409 410 psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); 411 412 e0->data.F32[i] = pol.e0; 413 e1->data.F32[i] = pol.e1; 414 e2->data.F32[i] = pol.e2; 415 } 416 417 pmTrend2D *trend = psf->params->data[PM_PAR_E0]; 418 if (trend->mode == PM_TREND_MAP) { 419 float errorFloor = 0.0; 420 float errorTotal = 0.0; 421 float errorTotalMin = FLT_MAX; 422 int entryMin = -1; 423 424 // XXX do I need to store the initial value of psfTry->mask and reset? 425 426 // check the fit residuals and increase Nx,Ny until the error is minimized 427 for (int i = 1; i < 10; i++) { 428 429 if (!pmPSFParamTrend (psf, i, &errorFloor, &errorTotal, psfTry->mask, x, y, e0, e1, e2, dz)) { 430 break; 431 } 432 433 // we do not have a good model for the error distributions on e0, e1, e2. 434 // use the bright end scatter from the constant fit to set the scale for the 435 // versions with 2D variation 436 if (i == 1) { 437 for (int i = 0; i < psfTry->sources->n; i++) { 438 if (dz) { 439 dz->data.F32[i] = errorFloor; 440 } 441 } 442 } 443 444 // store the resulting errorTotal values and the scales, redo the best 445 if (errorTotal < errorTotalMin) { 446 errorTotalMin = errorTotal; 447 entryMin = i; 448 } 449 } 450 if (entryMin == -1) { 451 psAbort ("failed on ApResid Trend"); 452 } 453 454 // XXX catch error condition 455 pmPSFParamTrend (psf, entryMin, &errorFloor, &errorTotal, psfTry->mask, x, y, e0, e1, e2, dz); 456 } else { 457 // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. 458 // This way, the parameters masked by one of the fits will be applied to the others 459 for (int i = 0; i < psf->psfTrendStats->clipIter; i++) { 460 461 // XXX we are using the same stats structure on each pass: do we need to re-init it? 462 pmTrend2DFit (psf->params->data[PM_PAR_E0], psfTry->mask, 0xff, x, y, e0, dz); 463 psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n); 464 pmTrend2DFit (psf->params->data[PM_PAR_E1], psfTry->mask, 0xff, x, y, e1, dz); 465 psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n); 466 pmTrend2DFit (psf->params->data[PM_PAR_E2], psfTry->mask, 0xff, x, y, e2, dz); 467 psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n); 468 } 469 } 470 471 // test dump of the psf parameters 472 if (psTraceGetLevel("psModules.objects") >= 4) { 473 FILE *f = fopen ("pol.dat", "w"); 474 fprintf (f, "# x y : e0obs e1obs e2obs : e0fit e1fit e2fit : mask\n"); 475 for (int i = 0; i < e0->n; i++) { 476 fprintf (f, "%f %f %f : %f %f %f : %f %f %f : %d\n", 477 x->data.F32[i], y->data.F32[i], flux->data.F32[i], 478 e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], 479 pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]), 480 pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]), 481 pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]), 482 psfTry->mask->data.U8[i]); 483 } 484 fclose (f); 485 } 486 487 psFree (e0); 488 psFree (e1); 489 psFree (e2); 363 } 364 365 if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) { 366 psFree(x); 367 psFree(y); 368 psFree(z); 369 return false; 370 } 490 371 491 372 // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY) 492 // XXXapply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters373 // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters 493 374 for (int i = 0; i < psf->params->n; i++) { 494 375 switch (i) { … … 508 389 for (int j = 0; j < psfTry->sources->n; j++) { 509 390 pmSource *source = psfTry->sources->data[j]; 510 if (source->modelEXT == NULL) 511 continue; 391 if (source->modelEXT == NULL) continue; 512 392 z->data.F32[j] = source->modelEXT->params->data.F32[i]; 513 393 } 514 394 515 pmTrend2D *trendE0 = psf->params->data[PM_PAR_E0];516 517 // re-create the pmTrend2D using the scales defined for E0518 395 psImageBinning *binning = psImageBinningAlloc(); 519 if (trendE0->mode == PM_TREND_MAP) { 520 binning->nXruff = trend->map->map->numCols; 521 binning->nYruff = trend->map->map->numRows; 522 } else { 523 binning->nXruff = trend->poly->nX; 524 binning->nYruff = trend->poly->nY; 525 } 396 binning->nXruff = psf->trendNx; 397 binning->nYruff = psf->trendNy; 526 398 binning->nXfine = psf->fieldNx; 527 399 binning->nYfine = psf->fieldNy; 528 psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); 529 psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); 530 400 401 if (psf->psfTrendMode == PM_TREND_MAP) { 402 psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); 403 psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); 404 } 405 406 // free existing trend, re-alloc 531 407 psFree (psf->params->data[i]); 532 psf->params->data[i] = pmTrend2DNoImageAlloc ( trendE0->mode, binning, psf->psfTrendStats);408 psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats); 533 409 psFree (binning); 534 410 … … 541 417 psFree(y); 542 418 psFree(z); 543 psFree(dz);544 419 return false; 545 420 } … … 552 427 for (int j = 0; j < psfTry->sources->n; j++) { 553 428 pmSource *source = psfTry->sources->data[j]; 554 if (source == NULL) 555 continue; 556 557 pmModel *model = source->modelEXT; 558 if (model == NULL) 559 continue; 560 561 pmModel *modelPSF = pmModelFromPSF (model, psf); 562 563 fprintf (f, "%f %f : ", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); 429 if (source == NULL) continue; 430 if (source->modelEXT == NULL) continue; 431 432 pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf); 433 434 fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]); 564 435 565 436 for (int i = 0; i < psf->params->n; i++) { 566 437 if (psf->params->data[i] == NULL) continue; 567 fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]);438 fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]); 568 439 } 569 fprintf (f, "%f %d\n", model->chisq, model->nIter);440 fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter); 570 441 571 442 psFree(modelPSF); … … 577 448 psFree (y); 578 449 psFree (z); 579 psFree (dz);580 psFree (flux);581 450 return true; 582 451 } 583 452 584 // XXX this function only works for MAP type -- extend to poly... 585 bool pmPSFParamTrend (pmPSF *psf, int scale, float *errorFloor, float *errorTotal, psVector *mask, psVector *x, psVector *y, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) { 453 454 bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) { 455 456 // we are doing a robust fit. after each pass, we drop points which are more deviant than 457 // three sigma. mask is currently updated for each pass. 458 459 // The shape parameters (SXX, SXY, SYY) are strongly coupled. We have to handle them very 460 // carefully. First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for 461 // each source and fit this set of parameters. These values are less tightly coupled, but 462 // are still inter-related. The fitted values do a good job of constraining the major axis 463 // and the position angle, but the minor axis is weakly measured. When we apply the PSF 464 // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape 465 // parameters, with the constraint that the minor axis must be greater than a minimum 466 // threshold. 467 468 // convert the measured source shape paramters to polarization terms 469 psVector *e0 = psVectorAlloc (sources->n, PS_TYPE_F32); 470 psVector *e1 = psVectorAlloc (sources->n, PS_TYPE_F32); 471 psVector *e2 = psVectorAlloc (sources->n, PS_TYPE_F32); 472 psVector *mag = psVectorAlloc (sources->n, PS_TYPE_F32); 473 for (int i = 0; i < sources->n; i++) { 474 pmSource *source = sources->data[i]; 475 if (source->modelEXT == NULL) continue; 476 // XXX I am relying on the fact that none of the masked sources 477 // have modelEXT set here. perhaps use the value of psfTry->mask instead? 478 479 psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); 480 481 e0->data.F32[i] = pol.e0; 482 e1->data.F32[i] = pol.e1; 483 e2->data.F32[i] = pol.e2; 484 485 float flux = source->modelEXT->params->data.F32[PM_PAR_I0]; 486 mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0; 487 } 488 489 if (psf->psfTrendMode == PM_TREND_MAP) { 490 float errorFloor = 0.0; 491 float errorTotal = 0.0; 492 float errorTotalMin = FLT_MAX; 493 int entryMin = -1; 494 495 psVector *dz = NULL; 496 psVector *mask = psVectorAlloc (sources->n, PS_TYPE_U8); 497 498 // check the fit residuals and increase Nx,Ny until the error is minimized 499 // pmPSFParamTrend increases the number along the longer of x or y 500 for (int i = 1; i < PS_MAX (psf->trendNx, psf->trendNy); i++) { 501 502 psVectorInit (mask, 0); 503 if (!pmPSFFitShapeParamsMap (psf, i, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) { 504 break; 505 } 506 507 // we do not have a good model for the error distributions on e0, e1, e2. 508 // use the bright end scatter from the constant fit to set the scale for the 509 // versions with 2D variation. potentially scale by poisson errors... 510 if (i == 1) { 511 // if (psf->poissonErrorsParams) do something else.. 512 dz = psVectorAlloc (sources->n, PS_TYPE_F32); 513 for (int i = 0; i < sources->n; i++) { 514 dz->data.F32[i] = errorFloor; 515 } 516 } 517 518 // store the resulting errorTotal values and the scales, redo the best 519 if (errorTotal < errorTotalMin) { 520 errorTotalMin = errorTotal; 521 entryMin = i; 522 } 523 } 524 if (entryMin == -1) { 525 psError (PS_ERR_UNKNOWN, true, "failed to find image map for shape params"); 526 return false; 527 } 528 529 // XXX supply the resulting mask values back to srcMask 530 psVectorInit (mask, 0); 531 if (!pmPSFFitShapeParamsMap (psf, entryMin, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) { 532 psAbort ("failed pmPSFFitShapeParamsMap on second pass?"); 533 } 534 535 pmTrend2D *trend = psf->params->data[PM_PAR_E0]; 536 psf->trendNx = trend->map->map->numCols; 537 psf->trendNy = trend->map->map->numRows; 538 539 psFree (mask); 540 psFree (dz); 541 } else { 542 543 // XXX iterate Nx, Ny based on scatter? 544 // XXX we force the x & y order to be the same 545 // modify the order to correspond to the actual number of matched stars: 546 int order = PS_MAX (psf->trendNx, psf->trendNy); 547 if ((sources->n < 15) && (order >= 3)) order = 2; 548 if ((sources->n < 11) && (order >= 2)) order = 1; 549 if ((sources->n < 8) && (order >= 1)) order = 0; 550 if ((sources->n < 3)) { 551 psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); 552 return false; 553 } 554 psf->trendNx = order; 555 psf->trendNy = order; 556 557 // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. 558 // This way, the parameters masked by one of the fits will be applied to the others 559 for (int i = 0; i < psf->psfTrendStats->clipIter; i++) { 560 561 // XXX we are using the same stats structure on each pass: do we need to re-init it? 562 bool status = true; 563 status &= pmTrend2DFit (psf->params->data[PM_PAR_E0], srcMask, 0xff, x, y, e0, NULL); 564 psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n); 565 status &= pmTrend2DFit (psf->params->data[PM_PAR_E1], srcMask, 0xff, x, y, e1, NULL); 566 psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n); 567 status &= pmTrend2DFit (psf->params->data[PM_PAR_E2], srcMask, 0xff, x, y, e2, NULL); 568 psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n); 569 570 if (!status) { 571 psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); 572 return false; 573 } 574 } 575 } 576 577 // test dump of the psf parameters 578 if (psTraceGetLevel("psModules.objects") >= 4) { 579 FILE *f = fopen ("pol.dat", "w"); 580 fprintf (f, "# x y : e0obs e1obs e2obs : e0fit e1fit e2fit : mask\n"); 581 for (int i = 0; i < e0->n; i++) { 582 fprintf (f, "%f %f : %f %f %f : %f %f %f : %d\n", 583 x->data.F32[i], y->data.F32[i], 584 e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], 585 pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]), 586 pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]), 587 pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]), 588 srcMask->data.U8[i]); 589 } 590 fclose (f); 591 } 592 593 psFree (e0); 594 psFree (e1); 595 psFree (e2); 596 psFree (mag); 597 return true; 598 } 599 600 // fit the shape variations as a psImageMap for the given scale factor 601 bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *errorFloor, float *errorTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) { 586 602 587 603 int Nx, Ny; … … 639 655 psf->psfTrendStats->clipIter = nIter; // restore default setting 640 656 641 # if (0)642 pmTrend2D *trend;643 trend = psf->params->data[PM_PAR_E0];644 testSaveImage (NULL, trend->map->map, "e0.fits");645 testSaveImage (NULL, trend->map->error, "e0d.fits");646 trend = psf->params->data[PM_PAR_E1];647 testSaveImage (NULL, trend->map->map, "e1.fits");648 testSaveImage (NULL, trend->map->error, "e1d.fits");649 trend = psf->params->data[PM_PAR_E2];650 testSaveImage (NULL, trend->map->map, "e2.fits");651 testSaveImage (NULL, trend->map->error, "e2d.fits");652 # endif653 654 657 // construct the fitted values and the residuals 655 658 psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], x, y); … … 700 703 // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin 701 704 int nGroup = PS_MAX (3*Nx*Ny, 10); 702 pmPSF ErrorFloor (errorFloor, dz, e0res, e1res, e2res, mask, nGroup);705 pmPSFShapeParamsErrors (errorFloor, mag, e0res, e1res, e2res, mask, nGroup); 703 706 704 707 *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum); … … 718 721 } 719 722 720 bool pmPSFErrorFloor (float *errorFloor, psVector *dMag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup) { 723 // calculate the minimum scatter of the parameters 724 bool pmPSFShapeParamsErrors (float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup) { 721 725 722 726 psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); 723 727 724 728 // measure the trend in bins with 10 values each; if < 10 total, use them all 725 int nBin = PS_MAX ( dMag->n / nGroup, 1);729 int nBin = PS_MAX (mag->n / nGroup, 1); 726 730 727 731 // output vectors for ApResid trend 728 732 psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32); 729 733 730 // use dMag to group the dMag and dap vectors731 psVector *index = psVectorSortIndex (NULL, dMag);732 733 // subset vectors for dMag and dap values within the given range734 // use mag to group parameters in sequence 735 psVector *index = psVectorSortIndex (NULL, mag); 736 737 // subset vectors for mag and dap values within the given range 734 738 psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); 735 739 psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); … … 740 744 for (int i = 0; i < dSo->n; i++) { 741 745 int j; 742 for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) {746 for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) { 743 747 int N = index->data.U32[n]; 744 748 dE0subset->data.F32[j] = e0res->data.F32[N];
Note:
See TracChangeset
for help on using the changeset viewer.
