Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 5592) +++ /trunk/psphot/Makefile (revision 5593) @@ -36,6 +36,4 @@ $(SRC)/psLibUtils.$(ARCH).o \ $(SRC)/psLine.$(ARCH).o \ -$(SRC)/psMinimize.$(ARCH).o \ -$(SRC)/psPolynomials.$(ARCH).o \ $(SRC)/psEllipse.$(ARCH).o \ $(SRC)/psModulesUtils.$(ARCH).o \ @@ -55,6 +53,4 @@ $(SRC)/pmModelGroup.$(ARCH).o: $(MODELS) - -# deprecated FITSOURCE = \ Index: /trunk/psphot/doc/notes.txt =================================================================== --- /trunk/psphot/doc/notes.txt (revision 5592) +++ /trunk/psphot/doc/notes.txt (revision 5593) @@ -1,4 +1,10 @@ Notes on psphot + +2005.11.16 + + I have made some fixes to make psphot work with pslib rel8_0. I + need to merge psphot with the rel8_0 version of the psMinimize and + pmObjects code. 2005.09.06 Index: /trunk/psphot/src/models/pmModel_GAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_GAUSS.c (revision 5592) +++ /trunk/psphot/src/models/pmModel_GAUSS.c (revision 5593) @@ -124,5 +124,5 @@ for (int i = 4; i < 7; i++) { psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]); + out[i] = psPolynomial2DEval (poly, out[2], out[3]); } return(true); Index: /trunk/psphot/src/models/pmModel_PGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_PGAUSS.c (revision 5592) +++ /trunk/psphot/src/models/pmModel_PGAUSS.c (revision 5593) @@ -140,5 +140,5 @@ for (int i = 4; i < 7; i++) { psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]); + out[i] = psPolynomial2DEval (poly, out[2], out[3]); } return(true); Index: /trunk/psphot/src/models/pmModel_QGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_QGAUSS.c (revision 5592) +++ /trunk/psphot/src/models/pmModel_QGAUSS.c (revision 5593) @@ -170,5 +170,5 @@ for (int i = 4; i < 8; i++) { psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]); + out[i] = psPolynomial2DEval (poly, out[2], out[3]); } return(true); Index: /trunk/psphot/src/models/pmModel_SGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_SGAUSS.c (revision 5592) +++ /trunk/psphot/src/models/pmModel_SGAUSS.c (revision 5593) @@ -300,5 +300,5 @@ for (int i = 4; i < 9; i++) { psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]); + out[i] = psPolynomial2DEval (poly, out[2], out[3]); } return(true); Index: /trunk/psphot/src/models/pmModel_TGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_TGAUSS.c (revision 5592) +++ /trunk/psphot/src/models/pmModel_TGAUSS.c (revision 5593) @@ -173,5 +173,5 @@ for (int i = 4; i < 8; i++) { psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]); + out[i] = psPolynomial2DEval (poly, out[2], out[3]); } return(true); Index: /trunk/psphot/src/pmObjects_EAM.c =================================================================== --- /trunk/psphot/src/pmObjects_EAM.c (revision 5592) +++ /trunk/psphot/src/pmObjects_EAM.c (revision 5593) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-17 19:04:33 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-25 23:13:43 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -999,5 +999,5 @@ int Ncr = 0; int Nsatstar = 0; - // psRegion allArray; + psRegion allArray = psRegionSet (0, 0, 0, 0); psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); @@ -1031,8 +1031,6 @@ // XXX EAM : can we use the value of SATURATE if mask is NULL? - // XXX EAM *** serious errors in psImageCountPixelMask: inconsistent with psRegion and subimage - // allArray = psRegionSet (1, tmpSrc->mask->numCols, 1, tmpSrc->mask->numRows); - // int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED); - int Nsatpix = psImageCountPixelMask_EAM (tmpSrc->mask, PSPHOT_MASK_SATURATED); + // allArray = psRegionSet (0, 0, 0, 0); + int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED); // saturated star (size consistent with PSF or larger) @@ -1648,5 +1646,5 @@ Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So)); - // psMinimizeLMChi2_EAM takes wt = 1/dY^2 + // psMinimizeLMChi2 takes wt = 1/dY^2 if (source->weight->data.F32[i][j] == 0) { yErr->data.F32[tmpCnt] = 0.0; @@ -1701,10 +1699,10 @@ psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); - fitStatus = psMinimizeLMChi2_EAM(myMin, covar, params, paramMask, x, y, yErr, modelFunc); + + fitStatus = psMinimizeLMChi2 (myMin, covar, params, paramMask, x, y, yErr, modelFunc); for (int i = 0; i < dparams->n; i++) { if ((paramMask != NULL) && paramMask->data.U8[i]) continue; dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); } - // XXX EAM: we need to do something (give an error?) if rc is false // XXX EAM: psMinimizeLMChi2 does not check convergence @@ -1724,5 +1722,5 @@ if (paramMask != NULL) { psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); - psMinimizeGaussNewtonDelta_EAM (delta, params, NULL, x, y, yErr, modelFunc); + psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc); for (int i = 0; i < dparams->n; i++) { if (!paramMask->data.U8[i]) continue; Index: /trunk/psphot/src/pmPSF.c =================================================================== --- /trunk/psphot/src/pmPSF.c (revision 5592) +++ /trunk/psphot/src/pmPSF.c (revision 5593) @@ -38,10 +38,9 @@ psf->params = psArrayAlloc (Nparams - 4); for (int i = 0; i < psf->params->n; i++) { - // XXX EAM : make this a user-defined value? - // XXX EAM : future version (0.7.0?) psf->params->data[i] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1); + // XXX EAM : make polynomial order user-specified? if (1) { - psf->params->data[i] = Polynomial2DAlloc_EAM(PS_POLYNOMIAL_ORD, 1, 1); + psf->params->data[i] = psPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD); } else { - psf->params->data[i] = Polynomial2DAlloc_EAM(PS_POLYNOMIAL_ORD, 0, 0); + psf->params->data[i] = psPolynomial2DAlloc(0, 0, PS_POLYNOMIAL_ORD); } } @@ -93,5 +92,5 @@ // XXX EAM : this is the expected API (cycle 7? cycle 8?) - psf->params->data[i] = VectorClipFitPolynomial2D_EAM (psf->params->data[i], stats, mask, 0xff, z, dz, x, y); + psf->params->data[i] = psVectorClipFitPolynomial2D (psf->params->data[i], stats, mask, 0xff, z, dz, x, y); // XXX EAM : drop this when above is implemented... Index: /trunk/psphot/src/pmPSFtry.c =================================================================== --- /trunk/psphot/src/pmPSFtry.c (revision 5592) +++ /trunk/psphot/src/pmPSFtry.c (revision 5593) @@ -265,17 +265,10 @@ } - // linear fit to rfBin, daBin - psPolynomial1D *poly = Polynomial1DAlloc_EAM (PS_POLYNOMIAL_ORD, 1); - - // XXX EAM : this is the intended API (cycle 7? cycle 8?) + // linear clipped fit to rfBin, daBin + psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD); psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - poly = VectorClipFitPolynomial1D_EAM (poly, fitstat, maskB, 1, daBin, NULL, rfBin); - - // poly = VectorFitPolynomial1D_EAM (poly, maskB, 1, daBin, NULL, rfBin); - - // XXX EAM : replace this when the above version is implemented - // poly = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL); - - psVector *daBinFit = Polynomial1DEvalVector_EAM (poly, rfBin); + poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin); + + psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin); psVector *daResid = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit); @@ -326,21 +319,8 @@ stats->clipIter = 3; - // linear fit to rfBin, daBin - psPolynomial1D *poly = Polynomial1DAlloc_EAM (PS_POLYNOMIAL_ORD, 1); - - // XXX EAM : this is the intended API (cycle 7? cycle 8?) - poly = VectorClipFitPolynomial1D_EAM (poly, stats, try->mask, PSFTRY_MASK_ALL, try->metric, NULL, rflux); - + // linear clipped fit to rfBin, daBin + psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD); + poly = psVectorClipFitPolynomial1D (poly, stats, try->mask, PSFTRY_MASK_ALL, try->metric, NULL, rflux); fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); - - // these stats should come back from the fit above... - // psVector *daFit = Polynomial1DEvalVector_EAM (poly, rflux); - // psVector *daResid = (psVector *) psBinaryOp (NULL, (void *) try->metric, "-", (void *) daFit); - - // stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV); - // stats->clipIter = 3; - // stats->clipSigma = 3; - - // stats = psVectorStats (stats, daResid, NULL, maskB, 1); try->psf->ApResid = poly->coeff[0]; Index: /trunk/psphot/src/pmPeaksSigmaLimit.c =================================================================== --- /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 5592) +++ /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 5593) @@ -16,5 +16,5 @@ psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32); - psImageSmooth_EAM (smooth, SIGMA, NSIGMA); + psImageSmooth (smooth, SIGMA, NSIGMA); psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot")); Index: /trunk/psphot/src/psLibUtils.c =================================================================== --- /trunk/psphot/src/psLibUtils.c (revision 5592) +++ /trunk/psphot/src/psLibUtils.c (revision 5593) @@ -1,95 +1,9 @@ # include // for strncasecmp +# include // for write # include # include "psLibUtils.h" -// XXX EAM : this is NOT included in the C99 headers ?? -FILE *fdopen(int fildes, const char *mode); - -// XXX EAM : these utility functions should be added back into PSLib - -# if (0) // now in psLib (v8) -static psHash *timers = NULL; - -// start/restart a named timer -bool psTimerStart (char *name) { - - psTime *start; - - if (timers == NULL) timers = psHashAlloc (16); - - start = psTimeGetNow (PS_TIME_UTC); - psHashAdd (timers, name, start); - psFree (start); - return (TRUE); -} - -bool psTimerClear (char *name) { - - bool status; - - if (name == NULL) return false; - - status = psHashRemove (timers, name); - return (status); -} - -// get current elapsed time on named timer -psF64 psTimerMark (char *name) { - - psTime *start; - psTime *mark; - psF64 delta; - - if (timers == NULL) return (0); - - start = psHashLookup (timers, name); - if (start == NULL) return (0); - - mark = psTimeGetNow (PS_TIME_UTC); - delta = psTimeDelta (mark, start); - psFree (mark); - // psFree (start); -- XXX is psHashLookup not psMemCopying? - - return (delta); -} - -void psTimerFree () { - - psFree (timers); - p_psTimeFinalize(); - return; -} -# endif - -# if (0) -// find the location of the specified argument -int psArgumentGet (int argc, char **argv, char *arg) { - - int i; - - for (i = 0; i < argc; i++) { - if (!strcmp(argv[i], arg)) - return (i); - } - - return ((int)NULL); -} - -// remove the specified argument (by location) -int psArgumentRemove (int N, int *argc, char **argv) { - - int i; - - if ((N != (int)NULL) && (N != 0)) { - (*argc)--; - for (i = N; i < *argc; i++) { - argv[i] = argv[i+1]; - } - } - - return (N); -} -# endif - +// XXX EAM : These two functions (psLogArguments and psTraceArguments) have been +// rolled together into a single function in psLib // we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae) // 2 = default, -v = 3, -vv = 4, -vvv = 5 @@ -151,326 +65,18 @@ } -# if (0) -// alternate implementation of this function from pmObjects.c -// now defined in psLib SDRS as psImageRow -psVector *psGetRowVectorFromImage(psImage *image, - psU32 row) -{ - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - - psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); - memcpy (tmpVector->data.F32, image->data.F32[row], image->numCols*sizeof(psF32)); - return(tmpVector); -} -# endif - -// XXX EAM : this is now in psLib -void psImageSmooth_EAM (psImage *image, float sigma, float Nsigma) { - - int Nx, Ny, Npixel, Nrange; - float factor, g, s; - psVector *temp; - - // relevant terms - Nrange = sigma*Nsigma + 0.5; - Npixel = 2*Nrange + 1; - factor = -0.5/(sigma*sigma); - - Nx = image->numCols; - Ny = image->numRows; - - // generate gaussian - psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_F32); - for (int i = -Nrange; i < Nrange + 1; i++) { - gaussnorm->data.F32[i+Nrange] = exp (factor*i*i); - } - psF32 *gauss = &gaussnorm->data.F32[Nrange]; - - // smooth in X direction - temp = psVectorAlloc (Nx, PS_TYPE_F32); - for (int j = 0; j < Ny; j++) { - psF32 *vi = image->data.F32[j]; - psF32 *vo = temp->data.F32; - for (int i = 0; i < Nx; i++) { - g = s = 0; - for (int n = -Nrange; n < Nrange + 1; n++) { - if (i+n < 0) continue; - if (i+n >= Nx) continue; - s += gauss[n]*vi[i+n]; - g += gauss[n]; - } - vo[i] = s / g; - } - memcpy (image->data.F32[j], temp->data.F32, Nx*sizeof(psF32)); - } - psFree (temp); - - // smooth in Y direction - temp = psVectorAlloc (image->numRows, PS_TYPE_F32); - for (int i = 0; i < Nx; i++) { - psF32 *vo = temp->data.F32; - psF32 **vi = image->data.F32; - for (int j = 0; j < Ny; j++) { - g = s = 0; - for (int n = -Nrange; n < Nrange + 1; n++) { - if (j+n < 0) continue; - if (j+n >= Ny) continue; - s += gauss[n]*vi[j+n][i]; - g += gauss[n]; - } - vo[j] = s / g; - } - // replace temp in image - for (int j = 0; j < Ny; j++) { - vi[j][i] = vo[j]; - } - } - psFree (temp); - psFree (gaussnorm); -} - -bool psImageInit (psImage *image,...) { - - va_list argp; - psU8 vU8; - psF32 vF32; - psF64 vF64; - - if (image == NULL) return (false); - - va_start (argp, image); - - switch (image->type.type) { - case PS_TYPE_U8: - vU8 = va_arg (argp, psU32); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.U8[iy][ix] = vU8; - } - } - break; - - case PS_TYPE_F32: - vF32 = va_arg (argp, psF64); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.F32[iy][ix] = vF32; - } - } - return (true); - - case PS_TYPE_F64: - vF64 = va_arg (argp, psF64); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.F64[iy][ix] = vF64; - } - } - return (true); - - default: - psAbort ("psphot.psUtils", "datatype %d not defined in psImageInit\n", image->type); - return (false); - } - return (false); -} - -// now in psLib (v8) -// count number of pixels with given mask value -// XXX EAM : version in psLib is broken -int psImageCountPixelMask_EAM (psImage *mask, psU8 value) -{ - int Npixels = 0; - - for (int i = 0; i < mask->numRows; i++) { - for (int j = 0; j < mask->numCols; j++) { - if (mask->data.U8[i][j] & value) { - Npixels ++; - } - } - } - return (Npixels); -} - -// define a square region centered on the given coordinate -// XXX EAM : this is now in psLib -# if (0) -psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) { - psRegion *region; - region = psRegionAlloc (x - radius, x + radius + 1, - y - radius, y + radius + 1); - return (region); -} -# endif - -// set actual region based on image parameters: -// compensate for negative upper limits -// XXX this is inconsistent: the coordindates should always be in the parent -// frame, which means the negative values should subtract from Nx,Ny of -// the parent, not the child. but, we don't carry the dimensions of the -// parent in the psImage container. for now, use the child Nx,Ny -// force range to be on this subimage -// XXX EAM : this needs to be changed to use psRegion rather than psRegion* -// XXX EAM : this is now in psLib -# if (0) -psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in) { - - // x0,y0, x1,y1 are in *parent* units - - if (out == NULL) { - out = psRegionAlloc(in->x0, in->x1, in->y0, in->y1); - } else { - *out = *in; - } - // XXX these are probably wrong (see above) - if (out->x1 <= 0) { - out->x1 = image->col0 + image->numCols + out->x1; - } - if (out->y1 <= 0) { - out->y1 = image->row0 + image->numRows + out->y1; - } - - // force the lower-limits to be on the child - out->x0 = PS_MAX(image->col0, out->x0); - out->y0 = PS_MAX(image->row0, out->y0); - - // force the upper-limits to be on the child - out->x1 = PS_MIN(image->col0 + image->numCols, out->x1); - out->y1 = PS_MIN(image->row0 + image->numRows, out->y1); - return (out); -} -# endif - -// mask the area contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) -void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) { - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (ix + image->col0 < region->x0) continue; - if (ix + image->col0 >= region->x1) continue; - if (iy + image->row0 < region->y0) continue; - if (iy + image->row0 >= region->y1) continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -// mask the area not contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) -void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) { - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (ix + image->col0 < region->x0) goto maskit; - if (ix + image->col0 >= region->x1) goto maskit; - if (iy + image->row0 < region->y0) goto maskit; - if (iy + image->row0 >= region->y1) goto maskit; - continue; - maskit: - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -// mask the area contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) -void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) { - - double dx, dy, r2, R2; - - R2 = PS_SQR(radius); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - dx = ix + image->col0 - x; - dy = iy + image->row0 - y; - r2 = PS_SQR(dx) + PS_SQR(dy); - if (r2 > R2) continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -// mask the area contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) -void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) { - - double dx, dy, r2, R2; - - R2 = PS_SQR(radius); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - dx = ix + image->col0 - x; - dy = iy + image->row0 - y; - r2 = PS_SQR(dx) + PS_SQR(dy); - if (r2 < R2) continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -# if (0) -psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) { - - int nBin = (upper - lower) / delta; - - psVector *out = psVectorAlloc (nBin, type); - - for (int i = 0; i < nBin; i++) { - out->data.F64[i] = lower + i * delta; - } - - return (out); -} -# endif - // XXX EAM a utility function bool p_psVectorPrintRow (int fd, psVector *a, char *name) { - FILE *f; - f = fdopen(fd, "a+"); - fprintf (f, "vector: %s\n", name); + char line[1024]; + + sprintf (line, "vector: %s\n", name); + write (fd, line, strlen(line)); for (int i = 0; i < a[0].n; i++) { - fprintf (f, "%f ", p_psVectorGetElementF64(a, i)); + sprintf (line, "%f ", p_psVectorGetElementF64(a, i)); + write (fd, line, strlen(line)); } - fprintf (f, "\n"); - fclose(f); + write (fd, "\n", 1); return (true); } -// XXX EAM is the use of fdopen here a good way to do this? Index: unk/psphot/src/psMinimize.c =================================================================== --- /trunk/psphot/src/psMinimize.c (revision 5592) +++ (revision ) @@ -1,302 +1,0 @@ -# include -# include "psLibUtils.h" - -// XXX EAM : this file contains my alternate implementations of psMinimizeLMChi2 -// XXX EAM this implementation of MinLM includes limits on params & dparams - -psBool p_psMinLM_GuessABP( - psImage *Alpha, - psVector *Beta, - psVector *Params, - const psImage *alpha, - const psVector *beta, - const psVector *params, - const psVector *paramMask, - const psVector *beta_lim, - const psVector *params_min, - const psVector *params_max, - psF64 lambda); - -psBool psMinimizeLMChi2_EAM(psMinimization *min, - psImage *covar, - psVector *params, - const psVector *paramMask, - const psArray *x, - const psVector *y, - const psVector *yWt, - psMinimizeLMChi2Func func) -{ - // PS_PTR_CHECK_NULL(min, NULL); - // PS_VECTOR_CHECK_NULL(params, NULL); - // PS_VECTOR_CHECK_EMPTY(params, NULL); - // PS_PTR_CHECK_NULL(x, NULL); - // PS_VECTOR_CHECK_NULL(y, NULL); - // PS_VECTOR_CHECK_EMPTY(y, NULL); - // PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL); - // PS_PTR_CHECK_NULL(func, NULL); - - // this function has test and current values for several things - // the current best value is in lower case - // the next guess value is in upper case - - // allocate internal arrays (current vs Guess) - psImage *alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); - psImage *Alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); - psVector *beta = psVectorAlloc (params->n, PS_TYPE_F64); - psVector *Beta = psVectorAlloc (params->n, PS_TYPE_F64); - psVector *Params = psVectorAlloc (params->n, PS_TYPE_F32); - psVector *dy = NULL; - psF64 Chisq = 0.0; - psF64 lambda = 0.001; - - psVector *beta_lim = NULL; - psVector *param_min = NULL; - psVector *param_max = NULL; - - // if we are provided a covar image, we expect to find these three vectors in first three rows - if (covar != NULL) { - beta_lim = psVectorAlloc (params->n, PS_TYPE_F32); - param_min = psVectorAlloc (params->n, PS_TYPE_F32); - param_max = psVectorAlloc (params->n, PS_TYPE_F32); - for (int i = 0; i < params->n; i++) { - beta_lim->data.F32[i] = covar->data.F64[0][i]; - param_min->data.F32[i] = covar->data.F64[1][i]; - param_max->data.F32[i] = covar->data.F64[2][i]; - } - psImageRecycle (covar, params->n, params->n, PS_TYPE_F64); - } - - - // why is this needed here??? the initial guess on params is provided by the user - Params = psVectorCopy (Params, params, PS_TYPE_F32); - - // the user provides the error or NULL. we need to convert - // to appropriate weights - if (yWt != NULL) { - dy = (psVector *) yWt; - } else { - dy = psVectorAlloc (y->n, PS_TYPE_F32); - psVectorInit (dy, 1.0); -// for (int i = 0; i < dy->n; i++) { -// dy->data.F32[i] = 1.0; -// } - } - - // calculate initial alpha and beta, set chisq (min->value) - min->value = p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func); - if (isnan(min->value)) { - min->iter = min->maxIter; - return (false); - } -# ifndef PS_NO_TRACE - // dump some useful info if trace is defined - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { - p_psImagePrint (psTraceGetDestination(), alpha, "alpha guess"); - p_psVectorPrint (psTraceGetDestination(), beta, "beta guess"); - p_psVectorPrint (psTraceGetDestination(), params, "params guess"); - } - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) { - p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess"); - } -# endif /* PS_NO_TRACE */ - - // iterate until the tolerance is reached, or give up - while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) { - - // set a new guess for Alpha, Beta, Params - p_psMinLM_GuessABP_EAM (Alpha, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, lambda); - - // measure linear model prediction - psF64 dLinear = p_psMinLM_dLinear (Beta, beta, lambda); - -# ifndef PS_NO_TRACE - // dump some useful info if trace is defined - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { - p_psImagePrint (psTraceGetDestination(), Alpha, "alpha guess"); - p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess"); - p_psVectorPrint (psTraceGetDestination(), Params, "params guess"); - } - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) { - p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess"); - } -# endif /* PS_NO_TRACE */ - - // calculate Chisq for new guess, update Alpha & Beta - Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, paramMask, x, y, dy, func); - - // XXX EAM alternate convergence criterion: - // compare the delta (min->value - Chisq) with the - // expected delta from the linear model (dLinear) - // accept new guess (if improvement), or increase lambda - psF64 rho = (min->value - Chisq) / dLinear; - - psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho); -# ifndef PS_NO_TRACE - // dump some useful info if trace is defined - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { - p_psImagePrint (psTraceGetDestination(), Alpha, "alpha guess"); - p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess"); - p_psVectorPrint (psTraceGetDestination(), Params, "params guess"); - } -# endif /* PS_NO_TRACE */ - - /* if (Chisq < min->value) { */ - if (rho > 0.0) { - min->lastDelta = (min->value - Chisq) / (dy->n - params->n); - min->value = Chisq; - alpha = psImageCopy (alpha, Alpha, PS_TYPE_F64); - beta = psVectorCopy (beta, Beta, PS_TYPE_F64); - params = psVectorCopy (params, Params, PS_TYPE_F32); - lambda *= 0.1; - } else { - lambda *= 10.0; - } - min->iter ++; - } - psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter); - - // construct & return the covariance matrix (if requested) - if (covar != NULL) { - p_psMinLM_GuessABP_EAM (covar, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, 0.0); - } - - // free the internal temporary data - psFree (alpha); - psFree (Alpha); - psFree (beta); - psFree (Beta); - psFree (Params); - if (yWt == NULL) { - psFree (dy); - } - - if (min->iter == min->maxIter) { - return (false); - } - return (true); -} - -// XXX EAM : can we use static copies of LUv, LUm, A? -psBool p_psMinLM_GuessABP_EAM (psImage *Alpha, - psVector *Beta, - psVector *Params, - const psImage *alpha, - const psVector *beta, - const psVector *params, - const psVector *paramMask, - const psVector *beta_lim, - const psVector *params_min, - const psVector *params_max, - psF64 lambda) -{ - -# define USE_LU_DECOMP 1 -# if (USE_LU_DECOMP) - psVector *LUv = NULL; - psImage *LUm = NULL; - psImage *A = NULL; - psF32 det; - - // LU decomposition version - psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using LUD version\n"); - - // set new guess values (creates matrix A) - A = psImageCopy (NULL, alpha, PS_TYPE_F64); - for (int j = 0; j < params->n; j++) { - if ((paramMask != NULL) && (paramMask->data.U8[j])) continue; - A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); - } - - // solve A*beta = Beta (Alpha = 1/A) - // these operations do not modify the input values (creates LUm, LUv) - LUm = psMatrixLUD (NULL, &LUv, A); - Beta = psMatrixLUSolve (Beta, LUm, beta, LUv); - Alpha = psMatrixInvert (Alpha, A, &det); - -# else - // gauss-jordan version - psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using Gauss-J version"); - - // set new guess values (creates matrix A) - Beta = psVectorCopy (Beta, beta, PS_TYPE_F64); - Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64); - for (int j = 0; j < params->n; j++) { - if ((paramMask != NULL) && (paramMask->data.U8[j])) continue; - Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); - } - - psGaussJordan (Alpha, Beta); -# endif - - // apply Beta to get new Params values - for (int j = 0; j < params->n; j++) { - if ((paramMask != NULL) && (paramMask->data.U8[j])) continue; - // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; - // continue; - // compare Beta to beta limits - if (beta_lim != NULL) { - if (fabs(Beta->data.F64[j]) > fabs(beta_lim->data.F32[j])) { - Beta->data.F64[j] = (Beta->data.F64[j] > 0) ? fabs(beta_lim->data.F32[j]) : -fabs(beta_lim->data.F32[j]); - } - } - Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; - // compare new params to param limits - if (params_max != NULL) { - Params->data.F32[j] = PS_MIN (Params->data.F32[j], params_max->data.F32[j]); - } - if (params_min != NULL) { - Params->data.F32[j] = PS_MAX (Params->data.F32[j], params_min->data.F32[j]); - } - } - -# if (USE_LU_DECOMP) - psFree (A); - psFree (LUm); - psFree (LUv); -# endif - - return true; -} - -bool psMinimizeGaussNewtonDelta_EAM (psVector *delta, - const psVector *params, - const psVector *paramMask, - const psArray *x, - const psVector *y, - const psVector *yWt, - psMinimizeLMChi2Func func) -{ - - // allocate internal arrays (current vs Guess) - psImage *alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); - psImage *Alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); - psVector *beta = psVectorAlloc (params->n, PS_TYPE_F64); - psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64); - psVector *dy = NULL; - - // the user provides the error or NULL. we need to convert - // to appropriate weights - if (yWt != NULL) { - dy = (psVector *) yWt; - } else { - dy = psVectorAlloc (y->n, PS_TYPE_F32); - psVectorInit (dy, 1.0); -// for (int i = 0; i < dy->n; i++) { -// dy->data.F32[i] = 1.0; -// } - } - - p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func); - p_psMinLM_GuessABP (Alpha, delta, Params, alpha, beta, params, paramMask, NULL, NULL, NULL, 0.0); - - psFree (alpha); - psFree (Alpha); - psFree (beta); - psFree (Params); - if (yWt == NULL) { - psFree (dy); - } - return (true); -} - - Index: unk/psphot/src/psPolynomials.c =================================================================== --- /trunk/psphot/src/psPolynomials.c (revision 5592) +++ (revision ) @@ -1,610 +1,0 @@ -# include -# include "psLibUtils.h" - -// XXX EAM : this file defines alternate versions of the polynomial fitting -// functions. these differ from psLib v.0.5.0 in the definition -// of nOrder (here we use nOrder in the mathematical sense, not -// like IRAF. ie, a linear function is of nOrder 1) -// psLib v 0.7.0 should be in sync with out defintion of nOrder -// If so, we can drop these functions and use the psLib versions. - -// write out the terms of the given 1D polynomial -void Polynomial1DDump_EAM (psPolynomial1D *poly) { - - for (int i = 0; i < poly->nX + 1; i++) { - fprintf (stderr, "x^%d : %g +/- %g\n", i, poly->coeff[i], poly->coeffErr[i]); - } -} - -static void Polynomial1DFree_EAM (psPolynomial1D *poly) { - - if (poly == NULL) return; - psFree (poly->coeff); - psFree (poly->coeffErr); - psFree (poly->mask); - return; -} - -// XXX EAM : use Nterm = Norder + 1 definition -// XXX EAM : should we provide both order and nterms in struct? -psPolynomial1D* Polynomial1DAlloc_EAM(psPolynomialType type, - psS32 nOrder) -{ - // PS_INT_CHECK_NON_NEGATIVE(nOrder, NULL); - - psS32 i = 0; - psS32 nTerm = nOrder + 1; - psPolynomial1D* newPoly = NULL; - - newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); - newPoly->type = type; - newPoly->nX = nOrder; - newPoly->coeff = (psF64 *)psAlloc(nTerm * sizeof(psF64)); - newPoly->coeffErr = (psF64 *)psAlloc(nTerm * sizeof(psF64)); - newPoly->mask = (psU8 *)psAlloc(nTerm * sizeof(psU8)); - for (i = 0; i < nTerm; i++) { - newPoly->coeff[i] = 0.0; - newPoly->coeffErr[i] = 0.0; - newPoly->mask[i] = 0; - } - psMemSetDeallocator(newPoly, (psFreeFunc) Polynomial1DFree_EAM); - return(newPoly); -} - -psF64 Polynomial1DEval_EAM(const psPolynomial1D* myPoly, psF64 x) -{ - psS32 loop_x = 0; - psF64 polySum = 0.0; - psF64 xSum = 1.0; - - for (loop_x = 0; loop_x < myPoly->nX + 1; loop_x++) { - if (myPoly->mask[loop_x] == 0) { - polySum += xSum * myPoly->coeff[loop_x]; - } - xSum *= x; - } - - return(polySum); -} - -psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, - const psVector *x) -{ - // PS_POLY_CHECK_NULL(myPoly, NULL); - // PS_VECTOR_CHECK_NULL(x, NULL); - // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL); - - psVector *tmp; - - tmp = psVectorAlloc(x->n, PS_TYPE_F64); - for (psS32 i=0;in;i++) { - tmp->data.F64[i] = Polynomial1DEval_EAM(myPoly, x->data.F64[i]); - } - - return(tmp); -} - -// XXX EAM : my alternate BuildSums1D -static psVector *BuildSums1D_EAM(psVector* sums, - psF64 x, - psS32 nTerm) -{ - psS32 nSum = 0; - psF64 xSum = 0.0; - - nSum = 2*nTerm; - if (sums == NULL) { - sums = psVectorAlloc(nSum, PS_TYPE_F64); - } - if (nSum > sums->n) { - sums = psVectorRealloc(sums, nSum); - } - - xSum = 1.0; - for (int i = 0; i < nSum; i++) { - sums->data.F64[i] = xSum; - xSum *= x; - } - return (sums); -} - -// XXX EAM : test version of 1d fitting -psPolynomial1D* VectorFitPolynomial1D_EAM(psPolynomial1D* myPoly, - psVector *mask, - psMaskType maskValue, - const psVector *y, - const psVector *yErr, - const psVector *x) -{ - // I think this is 1 dimension down - psImage* A = NULL; - psVector* B = NULL; - psVector* xSums = NULL; - psS32 nTerm; - psF64 wt; - - psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4, - "---- VectorFitPolynomial1DOrd() begin ----\n"); - - // dump minutiae -# ifndef PS_NO_TRACE -# if (0) - // XXX EAM : change from FILE to fd breaks this code: - if (psTraceGetLevel (".psLib.dataManip.VectorFitPolynomial1DOrd") >= 5) { - FILE *f = psTraceGetDestination (); - fprintf (f, "VectorFitPolynomial1D()\n"); - for (int i = 0; i < x->n; i++) { - fprintf (f, "(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]); - } - } -# endif -# endif - - nTerm = myPoly->nX + 1; - A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); - B = psVectorAlloc(nTerm, PS_TYPE_F64); - - // Initialize data structures (why is this not a function!) - for (int i = 0; i < nTerm; i++) { - B->data.F64[i] = 0.0; - for (int j = 0; j < nTerm; j++) { - A->data.F64[i][j] = 0.0; - } - } - - // xSums look like: 1, x, x^2, ... x^(2n+1) - - // Build the B and A data structs. - for (int k = 0; k < x->n; k++) { - if ((mask != NULL) && (mask->data.U8[k] & maskValue)) continue; - xSums = BuildSums1D_EAM(xSums, x->data.F64[k], nTerm); - - if (yErr == NULL) { - wt = 1.0; - } else { - // this should probably by yErr^2 !! - // this should filter yErr == 0 values - wt = 1.0 / PS_SQR(yErr->data.F64[k]); - } - for (int i = 0; i < nTerm; i++) { - B->data.F64[i] += y->data.F64[k] * xSums->data.F64[i] * wt; - } - - // we could skip half of the array and assign at the end - // we must handle masked orders - for (int i = 0; i < nTerm; i++) { - for (int j = 0; j < nTerm; j++) { - A->data.F64[i][j] += xSums->data.F64[i + j] * wt; - } - } - } - - // GaussJordan version - if (0) { - // does the solution in place - psGaussJordan (A, B); - - // the first nTerm entries in B correspond directly to the desired - // polynomial coefficients. this is only true for the 1D case - for (int k = 0; k < nTerm; k++) { - myPoly->coeff[k] = B->data.F64[k]; - } - } - else - // LUD version of the fit - { - psImage *ALUD = NULL; - psVector* outPerm = NULL; - psVector* coeffs = NULL; - - ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); - ALUD = psMatrixLUD(ALUD, &outPerm, A); - coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - for (int k = 0; k < nTerm; k++) { - myPoly->coeff[k] = coeffs->data.F64[k]; - } - psFree (ALUD); - psFree (outPerm); - psFree (coeffs); - } - - psFree(A); - psFree(B); - psFree(xSums); - - psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4, - "---- VectorFitPolynomial1DOrd() begin ----\n"); - return (myPoly); -} - -// ********************** 2D polynomial functions ****************** - -static void Polynomial2DFree_EAM(psPolynomial2D *poly) { - - if (poly == NULL) return; - for (int i = 0; i < poly->nX + 1; i++) { - psFree (poly->coeff[i]); - psFree (poly->coeffErr[i]); - psFree (poly->mask[i]); - } - psFree (poly->coeff); - psFree (poly->coeffErr); - psFree (poly->mask); - return; -} - -// XXX EAM : use Nterm = Norder + 1 definition -// the user requests a polynomial of order Norder -psPolynomial2D* Polynomial2DAlloc_EAM(psPolynomialType type, - psS32 nXorder, - psS32 nYorder) -{ - // PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL); - // PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL); - - psS32 x = 0; - psS32 y = 0; - psS32 nXterm = nXorder + 1; - psS32 nYterm = nYorder + 1; - psPolynomial2D* newPoly = NULL; - - newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D)); - newPoly->type = type; - newPoly->nX = nXorder; - newPoly->nY = nYorder; - - newPoly->coeff = (psF64 **)psAlloc(nXterm * sizeof(psF64 *)); - newPoly->coeffErr = (psF64 **)psAlloc(nXterm * sizeof(psF64 *)); - newPoly->mask = (psMaskType **)psAlloc(nXterm * sizeof(psMaskType *)); - - for (x = 0; x < nXterm; x++) { - newPoly->coeff[x] = (psF64 *)psAlloc(nYterm * sizeof(psF64)); - newPoly->coeffErr[x] = (psF64 *)psAlloc(nYterm * sizeof(psF64)); - newPoly->mask[x] = (psU8 *)psAlloc(nYterm * sizeof(char)); - } - for (x = 0; x < nXterm; x++) { - for (y = 0; y < nYterm; y++) { - newPoly->coeff[x][y] = 0.0; - newPoly->coeffErr[x][y] = 0.0; - newPoly->mask[x][y] = 0; - } - } - psMemSetDeallocator(newPoly, (psFreeFunc) Polynomial2DFree_EAM); - return(newPoly); -} - -// write out the terms of the given 2D polynomial -void Polynomial2DDump_EAM (psPolynomial2D *poly) { - - for (int i = 0; i < poly->nX + 1; i++) { - for (int j = 0; j < poly->nY + 1; j++) { - fprintf (stderr, "x^%d y^%d : %g +/- %g\n", i, j, poly->coeff[i][j], poly->coeffErr[i][j]); - } - } -} - -psF64 Polynomial2DEval_EAM(const psPolynomial2D* myPoly, - psF64 x, - psF64 y) -{ - // PS_POLY_CHECK_NULL(myPoly, NAN); - - psS32 loop_x = 0; - psS32 loop_y = 0; - psF64 polySum = 0.0; - psF64 xSum = 1.0; - psF64 ySum = 1.0; - - // XXX EAM : nX is order, not nTerms - for (loop_x = 0; loop_x < myPoly->nX + 1; loop_x++) { - ySum = xSum; - // XXX EAM : nX is order, not nTerms - for (loop_y = 0; loop_y < myPoly->nY + 1; loop_y++) { - if (myPoly->mask[loop_x][loop_y] == 0) { - polySum += ySum * myPoly->coeff[loop_x][loop_y]; - } - ySum *= y; - } - xSum *= x; - } - - return(polySum); -} - -// XXX EAM : this version uses myPoly->nX as Norder, not Nterms -psVector *Polynomial2DEvalVector_EAM(const psPolynomial2D *myPoly, - const psVector *x, - const psVector *y) - -{ - // PS_POLY_CHECK_NULL(myPoly, NULL); - // PS_VECTOR_CHECK_NULL(x, NULL); - // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL); - // PS_VECTOR_CHECK_NULL(y, NULL); - // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL); - - psVector *tmp; - psS32 vecLen=x->n; - - // Determine the length of the output vector to by the minimum of the x,y vectors - if (y->n < vecLen) { - vecLen = y->n; - } - - // Create output vector to return - tmp = psVectorAlloc(vecLen, PS_TYPE_F32); - - // Evaluate the polynomial at the specified points - for (psS32 i=0; idata.F32[i] = Polynomial2DEval_EAM(myPoly, x->data.F32[i], y->data.F32[i]); - } - - // Return output vector - return(tmp); -} - -// XXX EAM : this version uses the F64 vectors -psVector *Polynomial2DEvalVectorD_EAM(const psPolynomial2D *myPoly, - const psVector *x, - const psVector *y) - -{ - // PS_POLY_CHECK_NULL(myPoly, NULL); - // PS_VECTOR_CHECK_NULL(x, NULL); - // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL); - // PS_VECTOR_CHECK_NULL(y, NULL); - // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL); - - psVector *tmp; - psS32 vecLen=x->n; - - // Determine the length of the output vector to by the minimum of the x,y vectors - if (y->n < vecLen) { - vecLen = y->n; - } - - // Create output vector to return - tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - - // Evaluate the polynomial at the specified points - for (psS32 i=0; idata.F64[i] = Polynomial2DEval_EAM(myPoly,x->data.F64[i],y->data.F64[i]); - } - - // Return output vector - return(tmp); -} - -// XXX EAM : BuildSums2D in analogy with BuildSums1D -static psImage *BuildSums2D_EAM(psImage* sums, - psF64 x, psF64 y, - psS32 nXterm, psS32 nYterm) -{ - psS32 nXsum = 0; - psS32 nYsum = 0; - psF64 xSum = 1.0; - psF64 ySum = 1.0; - - nXsum = 2*nXterm; - nYsum = 2*nYterm; - if (sums == NULL) { - sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64); - } - if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) { - psFree (sums); - sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64); - } - - ySum = 1.0; - for (int j = 0; j < nYsum; j++) { - xSum = ySum; - for (int i = 0; i < nXsum; i++) { - sums->data.F64[i][j] = xSum; - xSum *= x; - } - ySum *= y; - } - return (sums); -} - -// XXX EAM : test version of 2d fitting -psPolynomial2D* VectorFitPolynomial2D_EAM(psPolynomial2D* myPoly, - psVector* mask, - psMaskType maskValue, - const psVector* z, - const psVector* zErr, - const psVector* x, - const psVector* y) -{ - // I think this is 1 dimension down - psImage* A = NULL; - psVector* B = NULL; - psImage* Sums = NULL; - psF64 wt; - psS32 nTerm, nXterm, nYterm; - - nXterm = myPoly->nX + 1; - nYterm = myPoly->nY + 1; - nTerm = nXterm * nYterm; - - A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); - B = psVectorAlloc(nTerm, PS_TYPE_F64); - - // Initialize data structures (why is this not a function!) - for (int i = 0; i < nTerm; i++) { - B->data.F64[i] = 0.0; - for (int j = 0; j < nTerm; j++) { - A->data.F64[i][j] = 0.0; - } - } - - // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1) - - // Build the B and A data structs. - for (int k = 0; k < x->n; k++) { - if ((mask != NULL) && (mask->data.U8[k] & maskValue)) continue; - Sums = BuildSums2D_EAM(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm); - - if (zErr == NULL) { - wt = 1.0; - } else { - // this should probably by zErr^2 !! - // this should filter zErr == 0 values - wt = 1.0 / zErr->data.F64[k]; - } - - // we could skip half of the array and assign at the end - // we must handle masked orders - for (int n = 0; n < nXterm; n++) { - for (int m = 0; m < nYterm; m++) { - B->data.F64[n+m*nXterm] += z->data.F64[k] * Sums->data.F64[n][m] * wt; - } - } - - for (int i = 0; i < nXterm; i++) { - for (int j = 0; j < nYterm; j++) { - for (int n = 0; n < nXterm; n++) { - for (int m = 0; m < nYterm; m++) { - A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt; - } - } - } - } - } - - // does the solution in place - psGaussJordan (A, B); - - // XXX: How do we know if these routines were successful? - // ALUD = psMatrixLUD(ALUD, &outPerm, A); - // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - - for (int n = 0; n < nXterm; n++) { - for (int m = 0; m < nYterm; m++) { - myPoly->coeff[n][m] = B->data.F64[n+m*nXterm]; - } - } - - psFree(A); - psFree(B); - psFree(Sums); - - psTrace(".psLib.dataManip.VectorFitPolynomial2DOrd", 4, - "---- VectorFitPolynomial2DOrd() begin ----\n"); - return (myPoly); -} - -// XXX EAM : be careful here with F32 vs F64 vectors -psPolynomial1D* VectorClipFitPolynomial1D_EAM(psPolynomial1D* poly, - psStats *stats, - psVector* mask, - psMaskType maskValue, - const psVector* z, - const psVector* dz, - const psVector* x) -{ - psVector *zFit = NULL; - psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64); - float minClipSigma = stats->min; - float maxClipSigma = stats->max; - float minClipValue; - float maxClipValue; - - // XXX EAM : use SAMPLE_MEDIAN and SAMPLE_STDEV for stats: - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - - for (int N = 0; N < stats->clipIter; N++) { - int Nkeep = 0; - - poly = VectorFitPolynomial1D_EAM (poly, mask, maskValue, z, dz, x); - zFit = Polynomial1DEvalVector_EAM (poly, x); - zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit); - - stats = psVectorStats (stats, zResid, NULL, mask, maskValue); - minClipValue = -1*fabs(minClipSigma)*stats->sampleStdev; - maxClipValue = +1*fabs(maxClipSigma)*stats->sampleStdev; - - psTrace (".psphot.VectorClipFit", 4, "residual stats for robust fit: %g +/- %g\n", - stats->sampleMedian, stats->sampleStdev); - - // set mask if pts are not valid - // we are masking out any point which is out of range - // recovery is not allowed with this scheme - for (int i = 0; i < zResid->n; i++) { - if (mask->data.U8[i]) continue; - if (zResid->data.F64[i] - stats->sampleMedian > maxClipValue) { - mask->data.U8[i] |= 0x01; - continue; - } - if (zResid->data.F64[i] - stats->sampleMedian < minClipValue) { - mask->data.U8[i] |= 0x01; - continue; - } - Nkeep ++; - } - - psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n", - Nkeep, x->n); - - psFree (zFit); - } - psFree (zResid); - return (poly); -} - - -// XXX EAM : be careful here with F32 vs F64 vectors -psPolynomial2D* VectorClipFitPolynomial2D_EAM(psPolynomial2D* poly, - psStats *stats, - psVector* mask, - psMaskType maskValue, - const psVector* z, - const psVector* dz, - const psVector* x, - const psVector* y) -{ - psVector *zFit = NULL; - psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64); - - // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats: - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - - for (int N = 0; N < 3; N++) { - int Nkeep = 0; - - poly = VectorFitPolynomial2D_EAM (poly, mask, maskValue, z, dz, x, y); - zFit = Polynomial2DEvalVectorD_EAM (poly, x, y); - zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit); - - stats = psVectorStats (stats, zResid, NULL, mask, maskValue); - psTrace (".psphot.RobustFit", 4, "residual stats for robust fit: %g +/- %g\n", - stats->sampleMedian, stats->sampleStdev); - - // set mask if pts are not valid - // we are masking out any point which is out of range - // recovery is not allowed with this scheme - for (int i = 0; i < zResid->n; i++) { - if (mask->data.U8[i]) continue; - if (fabs(zResid->data.F64[i] - stats->sampleMedian) > 3*stats->sampleStdev) { - mask->data.U8[i] |= 0x01; - continue; - } - Nkeep ++; - } - - psTrace (".psphot.RobustFit", 4, "keeping %d of %d pts for fit\n", - Nkeep, x->n); - -// psTrace (".psphot.RobustFit", 4, "model pars: %f %f %f %f\n", -// poly->coeff[0][0], poly->coeff[1][0], -// poly->coeff[0][1], poly->coeff[1][1]); - - psTrace (".psphot.RobustFit", 4, "model pars: %f\n", poly->coeff[0][0]); - - psFree (zFit); - } - psFree (zResid); - return (poly); -} - -// XXX EAM : VectorFitPolynomial2D and Polynomial2DEvalVector require different types (F32 vs F64) Index: /trunk/psphot/src/psphotFitGalaxies.c =================================================================== --- /trunk/psphot/src/psphotFitGalaxies.c (revision 5592) +++ /trunk/psphot/src/psphotFitGalaxies.c (revision 5593) @@ -75,6 +75,4 @@ // fit FLT (not PSF) model (set/unset the pixel mask) psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); - // XXX EAM : 0.5.0 code used _EAM version - // status = pmSourceFitModel (source, model, false); status = pmSourceFitModel (source, model, false); psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);