Index: trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- trunk/psModules/src/objects/pmPSFtry.c (revision 15562) +++ trunk/psModules/src/objects/pmPSFtry.c (revision 15697) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.50 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-10 01:09:20 $ + * @version $Revision: 1.51 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-27 03:14:57 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -56,5 +56,5 @@ // allocate a pmPSFtry based on the desired sources and the model (identified by name) -pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options) +pmPSFtry *pmPSFtryAlloc (const psArray *sources, const pmPSFOptions *options) { pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry)); @@ -98,5 +98,5 @@ // generate a pmPSFtry with a copy of the test PSF sources -pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, pmPSFOptions *options, psMaskType maskVal, psMaskType mark) +pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psMaskType maskVal, psMaskType mark) { bool status; @@ -377,8 +377,8 @@ if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) { - psFree(x); - psFree(y); - psFree(z); - return false; + psFree(x); + psFree(y); + psFree(z); + return false; } @@ -406,19 +406,19 @@ } - psImageBinning *binning = psImageBinningAlloc(); - binning->nXruff = psf->trendNx; - binning->nYruff = psf->trendNy; - binning->nXfine = psf->fieldNx; - binning->nYfine = psf->fieldNy; - - if (psf->psfTrendMode == PM_TREND_MAP) { - psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); - psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); - } - - // free existing trend, re-alloc - psFree (psf->params->data[i]); - psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats); - psFree (binning); + psImageBinning *binning = psImageBinningAlloc(); + binning->nXruff = psf->trendNx; + binning->nYruff = psf->trendNy; + binning->nXfine = psf->fieldNx; + binning->nYfine = psf->fieldNy; + + if (psf->psfTrendMode == PM_TREND_MAP) { + psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); + } + + // free existing trend, re-alloc + psFree (psf->params->data[i]); + psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats); + psFree (binning); // fit the collection of measured parameters to the PSF 2D model @@ -487,6 +487,6 @@ pmSource *source = sources->data[i]; if (source->modelEXT == NULL) continue; - // XXX I am relying on the fact that none of the masked sources - // have modelEXT set here. perhaps use the value of psfTry->mask instead? + // XXX I am relying on the fact that none of the masked sources + // have modelEXT set here. perhaps use the value of psfTry->mask instead? psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); @@ -496,94 +496,94 @@ e2->data.F32[i] = pol.e2; - float flux = source->modelEXT->params->data.F32[PM_PAR_I0]; - mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0; + float flux = source->modelEXT->params->data.F32[PM_PAR_I0]; + mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0; } if (psf->psfTrendMode == PM_TREND_MAP) { - float errorFloor = 0.0; - float errorTotal = 0.0; - float errorTotalMin = FLT_MAX; - int entryMin = -1; - - psVector *dz = NULL; - psVector *mask = psVectorAlloc (sources->n, PS_TYPE_U8); - - // check the fit residuals and increase Nx,Ny until the error is minimized - // pmPSFParamTrend increases the number along the longer of x or y - for (int i = 1; i < PS_MAX (psf->trendNx, psf->trendNy); i++) { - - psVectorInit (mask, 0); - if (!pmPSFFitShapeParamsMap (psf, i, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) { - break; - } - - // we do not have a good model for the error distributions on e0, e1, e2. - // use the bright end scatter from the constant fit to set the scale for the - // versions with 2D variation. potentially scale by poisson errors... - if (i == 1) { - // if (psf->poissonErrorsParams) do something else.. - dz = psVectorAlloc (sources->n, PS_TYPE_F32); - for (int i = 0; i < sources->n; i++) { - dz->data.F32[i] = errorFloor; - } - } - - // store the resulting errorTotal values and the scales, redo the best - if (errorTotal < errorTotalMin) { - errorTotalMin = errorTotal; - entryMin = i; - } - } - if (entryMin == -1) { - psError (PS_ERR_UNKNOWN, true, "failed to find image map for shape params"); - return false; - } - - // XXX supply the resulting mask values back to srcMask - psVectorInit (mask, 0); - if (!pmPSFFitShapeParamsMap (psf, entryMin, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) { - psAbort ("failed pmPSFFitShapeParamsMap on second pass?"); - } - - pmTrend2D *trend = psf->params->data[PM_PAR_E0]; - psf->trendNx = trend->map->map->numCols; - psf->trendNy = trend->map->map->numRows; - - psFree (mask); - psFree (dz); + float errorFloor = 0.0; + float errorTotal = 0.0; + float errorTotalMin = FLT_MAX; + int entryMin = -1; + + psVector *dz = NULL; + psVector *mask = psVectorAlloc (sources->n, PS_TYPE_U8); + + // check the fit residuals and increase Nx,Ny until the error is minimized + // pmPSFParamTrend increases the number along the longer of x or y + for (int i = 1; i < PS_MAX (psf->trendNx, psf->trendNy); i++) { + + psVectorInit (mask, 0); + if (!pmPSFFitShapeParamsMap (psf, i, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) { + break; + } + + // we do not have a good model for the error distributions on e0, e1, e2. + // use the bright end scatter from the constant fit to set the scale for the + // versions with 2D variation. potentially scale by poisson errors... + if (i == 1) { + // if (psf->poissonErrorsParams) do something else.. + dz = psVectorAlloc (sources->n, PS_TYPE_F32); + for (int i = 0; i < sources->n; i++) { + dz->data.F32[i] = errorFloor; + } + } + + // store the resulting errorTotal values and the scales, redo the best + if (errorTotal < errorTotalMin) { + errorTotalMin = errorTotal; + entryMin = i; + } + } + if (entryMin == -1) { + psError (PS_ERR_UNKNOWN, true, "failed to find image map for shape params"); + return false; + } + + // XXX supply the resulting mask values back to srcMask + psVectorInit (mask, 0); + if (!pmPSFFitShapeParamsMap (psf, entryMin, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) { + psAbort ("failed pmPSFFitShapeParamsMap on second pass?"); + } + + pmTrend2D *trend = psf->params->data[PM_PAR_E0]; + psf->trendNx = trend->map->map->numCols; + psf->trendNy = trend->map->map->numRows; + + psFree (mask); + psFree (dz); } else { - // XXX iterate Nx, Ny based on scatter? - // XXX we force the x & y order to be the same - // modify the order to correspond to the actual number of matched stars: - int order = PS_MAX (psf->trendNx, psf->trendNy); - if ((sources->n < 15) && (order >= 3)) order = 2; - if ((sources->n < 11) && (order >= 2)) order = 1; - if ((sources->n < 8) && (order >= 1)) order = 0; - if ((sources->n < 3)) { - psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); - return false; - } - psf->trendNx = order; - psf->trendNy = order; - - // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. - // This way, the parameters masked by one of the fits will be applied to the others - for (int i = 0; i < psf->psfTrendStats->clipIter; i++) { - - // XXX we are using the same stats structure on each pass: do we need to re-init it? - bool status = true; - status &= pmTrend2DFit (psf->params->data[PM_PAR_E0], srcMask, 0xff, x, y, e0, NULL); - psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n); - status &= pmTrend2DFit (psf->params->data[PM_PAR_E1], srcMask, 0xff, x, y, e1, NULL); - psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n); - status &= pmTrend2DFit (psf->params->data[PM_PAR_E2], srcMask, 0xff, x, y, e2, NULL); - psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n); - - if (!status) { - psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); - return false; - } - } + // XXX iterate Nx, Ny based on scatter? + // XXX we force the x & y order to be the same + // modify the order to correspond to the actual number of matched stars: + int order = PS_MAX (psf->trendNx, psf->trendNy); + if ((sources->n < 15) && (order >= 3)) order = 2; + if ((sources->n < 11) && (order >= 2)) order = 1; + if ((sources->n < 8) && (order >= 1)) order = 0; + if ((sources->n < 3)) { + psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); + return false; + } + psf->trendNx = order; + psf->trendNy = order; + + // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. + // This way, the parameters masked by one of the fits will be applied to the others + for (int i = 0; i < psf->psfTrendStats->clipIter; i++) { + + // XXX we are using the same stats structure on each pass: do we need to re-init it? + bool status = true; + status &= pmTrend2DFit (psf->params->data[PM_PAR_E0], srcMask, 0xff, x, y, e0, NULL); + psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n); + status &= pmTrend2DFit (psf->params->data[PM_PAR_E1], srcMask, 0xff, x, y, e1, NULL); + psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n); + status &= pmTrend2DFit (psf->params->data[PM_PAR_E2], srcMask, 0xff, x, y, e2, NULL); + psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n); + + if (!status) { + psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); + return false; + } + } } @@ -591,13 +591,13 @@ if (psTraceGetLevel("psModules.objects") >= 4) { FILE *f = fopen ("pol.dat", "w"); - fprintf (f, "# x y : e0obs e1obs e2obs : e0fit e1fit e2fit : mask\n"); + fprintf (f, "# x y : e0obs e1obs e2obs : e0fit e1fit e2fit : mask\n"); for (int i = 0; i < e0->n; i++) { fprintf (f, "%f %f : %f %f %f : %f %f %f : %d\n", x->data.F32[i], y->data.F32[i], - e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], + e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]), pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]), pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]), - srcMask->data.U8[i]); + srcMask->data.U8[i]); } fclose (f); @@ -615,17 +615,17 @@ int Nx, Ny; - + // set the map scale to match the aspect ratio if (psf->fieldNx > psf->fieldNy) { - Nx = scale; - Ny = PS_MAX (1, Nx * (psf->fieldNy / psf->fieldNx)); + Nx = scale; + Ny = PS_MAX (1, Nx * (psf->fieldNy / psf->fieldNx)); } else { - Ny = scale; - Nx = PS_MAX (1, Ny * (psf->fieldNx / psf->fieldNy)); + Ny = scale; + Nx = PS_MAX (1, Ny * (psf->fieldNx / psf->fieldNy)); } // do we have enough sources for this fine of a grid? if (x->n < 3*Nx*Ny) { - return false; + return false; } @@ -658,5 +658,5 @@ // This way, the parameters masked by one of the fits will be applied to the others for (int i = 0; i < nIter; i++) { - // XXX we are using the same stats structure on each pass: do we need to re-init it? + // XXX we are using the same stats structure on each pass: do we need to re-init it? pmTrend2DFit (psf->params->data[PM_PAR_E0], mask, 0xff, x, y, e0obs, dz); psTrace ("psModules.objects", 5, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0obs->n); @@ -683,8 +683,8 @@ float mapErrorSum = 0.0; float mapError = 0.0; - for (int iy = 0; iy < trend->map->error->numRows; iy++) { - for (int ix = 0; ix < trend->map->error->numCols; ix++) { - mapError += PS_SQR (trend->map->error->data.F32[iy][ix]); - } + for (int iy = 0; iy < trend->map->error->numRows; iy++) { + for (int ix = 0; ix < trend->map->error->numCols; ix++) { + mapError += PS_SQR (trend->map->error->data.F32[iy][ix]); + } } psTrace ("psModules.objects", 4, "E0 error: %f\n", sqrt(mapError)); @@ -693,8 +693,8 @@ trend = psf->params->data[PM_PAR_E1]; mapError = 0.0; - for (int iy = 0; iy < trend->map->error->numRows; iy++) { - for (int ix = 0; ix < trend->map->error->numCols; ix++) { - mapError += PS_SQR (trend->map->error->data.F32[iy][ix]); - } + for (int iy = 0; iy < trend->map->error->numRows; iy++) { + for (int ix = 0; ix < trend->map->error->numCols; ix++) { + mapError += PS_SQR (trend->map->error->data.F32[iy][ix]); + } } psTrace ("psModules.objects", 4, "E1 error: %f\n", sqrt(mapError)); @@ -703,8 +703,8 @@ trend = psf->params->data[PM_PAR_E2]; mapError = 0.0; - for (int iy = 0; iy < trend->map->error->numRows; iy++) { - for (int ix = 0; ix < trend->map->error->numCols; ix++) { - mapError += PS_SQR (trend->map->error->data.F32[iy][ix]); - } + for (int iy = 0; iy < trend->map->error->numRows; iy++) { + for (int ix = 0; ix < trend->map->error->numCols; ix++) { + mapError += PS_SQR (trend->map->error->data.F32[iy][ix]); + } } psTrace ("psModules.objects", 4, "E2 error: %f\n", sqrt(mapError)); @@ -742,5 +742,5 @@ int nBin = PS_MAX (mag->n / nGroup, 1); - // output vectors for ApResid trend + // output vectors for ApResid trend psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32); @@ -756,33 +756,33 @@ int n = 0; for (int i = 0; i < dSo->n; i++) { - int j; - for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) { - int N = index->data.U32[n]; - dE0subset->data.F32[j] = e0res->data.F32[N]; - dE1subset->data.F32[j] = e1res->data.F32[N]; - dE2subset->data.F32[j] = e2res->data.F32[N]; - - mkSubset->data.U8[j] = mask->data.U8[N]; - } - dE0subset->n = j; - dE1subset->n = j; - dE2subset->n = j; - mkSubset->n = j; - - // calculate the root-mean-square of E0, E1, E2 - float dEsquare = 0.0; - psStatsInit (statsS); - psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff); - dEsquare += PS_SQR(statsS->robustStdev); - - psStatsInit (statsS); - psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff); - dEsquare += PS_SQR(statsS->robustStdev); - - psStatsInit (statsS); - psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff); - dEsquare += PS_SQR(statsS->robustStdev); - - dSo->data.F32[i] = sqrt(dEsquare); + int j; + for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) { + int N = index->data.U32[n]; + dE0subset->data.F32[j] = e0res->data.F32[N]; + dE1subset->data.F32[j] = e1res->data.F32[N]; + dE2subset->data.F32[j] = e2res->data.F32[N]; + + mkSubset->data.U8[j] = mask->data.U8[N]; + } + dE0subset->n = j; + dE1subset->n = j; + dE2subset->n = j; + mkSubset->n = j; + + // calculate the root-mean-square of E0, E1, E2 + float dEsquare = 0.0; + psStatsInit (statsS); + psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff); + dEsquare += PS_SQR(statsS->robustStdev); + + psStatsInit (statsS); + psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff); + dEsquare += PS_SQR(statsS->robustStdev); + + psStatsInit (statsS); + psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff); + dEsquare += PS_SQR(statsS->robustStdev); + + dSo->data.F32[i] = sqrt(dEsquare); } psFree (dE0subset); @@ -790,5 +790,5 @@ psFree (dE2subset); psFree (mkSubset); - + psStats *stats = psStatsAlloc (PS_STAT_MIN); psVectorStats (stats, dSo, NULL, NULL, 0); @@ -800,5 +800,5 @@ psFree (dSo); - + psFree (statsS);