Index: /branches/eam_branches/ipp-20250626/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20250626/psModules/src/objects/Makefile.am (revision 42954) +++ /branches/eam_branches/ipp-20250626/psModules/src/objects/Makefile.am (revision 42955) @@ -66,4 +66,6 @@ pmSourcePlotMoments.c \ pmSourcePlotApResid.c \ + pmSourceTrace.c \ + pmSourceTraceFitGauss.c \ pmSourceVisual.c \ pmResiduals.c \ @@ -125,4 +127,5 @@ pmSourcePlots.h \ pmSourceVisual.h \ + pmSourceTrace.h \ pmResiduals.h \ pmPSF.h \ Index: /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSource.h =================================================================== --- /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSource.h (revision 42954) +++ /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSource.h (revision 42955) @@ -27,4 +27,5 @@ PM_SOURCE_TYPE_STAR, ///< a good-quality star (subtracted model is PSF) PM_SOURCE_TYPE_EXTENDED, ///< an extended object (eg, galaxy) (subtracted model is EXT) + PM_SOURCE_TYPE_TRACE, ///< a 1D profile x Gaussian (other) PSF } pmSourceType; Index: /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSourceTrace.c =================================================================== --- /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSourceTrace.c (revision 42955) +++ /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSourceTrace.c (revision 42955) @@ -0,0 +1,262 @@ +/** @file pmSourceTrace.c + * @author EAM, IfA + * @date $Date: 2008-10-08 21:53:08 $ + * Copyright 2025 IfA, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmSourceTrace.h" + +/******************************************************************************/ +static void pmSourceTraceFree(pmSourceTrace *tmp) +{ + if (!tmp) return; + psFree (tmp->profile); + return; +} + +// define a source trace with a specified origin, position angle, sigma +// profile is undefined to start +pmSourceTrace *pmSourceTraceAlloc(psF32 x, psF32 y, psF32 posAngle, psF32 sigma) { + + psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__); + pmSourceTrace *tmp = (pmSourceTrace *) psAlloc(sizeof(pmSourceTrace)); + + tmp->profile = NULL; + + tmp->x = x; + tmp->y = y; + tmp->posAngle = posAngle; + tmp->sigma = sigma; + + psMemSetDeallocator(tmp, (psFreeFunc) pmSourceTraceFree); + + psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__); + return(tmp); +} + +# define dCOS(DEGREES) (cos(DEGREES * PS_RAD_DEG)) +# define dSIN(DEGREES) (sin(DEGREES * PS_RAD_DEG)) + +// insert source in the image +bool pmSourceTraceAdd (psImage *image, pmSourceTrace *source) { + + if (!image) return false; + if (!source) return false; + if (!source->profile) return false; + + // inefficient method: run through pixels in the output image + // and calculate model flux at that location + + float CS = dCOS(source->posAngle); + float SN = dSIN(source->posAngle); + float sigma = source->sigma; + + int Npix = source->profile->n; + int Nx = image->numCols; + int Ny = image->numRows; + + if (Npix < 2) return false; + + for (int iy = 0; iy < Ny; iy++) { + for (int ix = 0; ix < Nx; ix++) { + + // transform (x,y) to (p,q) where p is along the profile an q is cross profile + // (p,q) = (0,0) corresponds to src->(x,y) + + // source->x,y are in pixel coordinates, not index + // convert ix,iy to pixel coords + float dx = ix + 0.5 - source->x; + float dy = iy + 0.5 - source->y; + + // p is the coordinate along the profile. + float p = CS*dx + SN*dy; + float q = CS*dy - SN*dx; + + // allow extrapolation 1/2 pixel beyond the end of the profile + if (p < -0.5) continue; + if (p > Npix - 0.5) continue; + + // the profile vector represents + // interpolate the nearest pixels + + // 5.5 is center of pixel with index 5. If coordinate is 5.5 - 6.5 + // interpolate between the pixel with index 5 and index 6. + // int (5.51 - 0.5) = int (5.01) = 5 + // int (6.49 - 0.5) = int (5.99) = 5 + // int (5.51 + 0.5) = int (6.01) = 6 + // int (6.49 - 0.5) = int (6.99) = 6 + // lower pixel = (int) (p - 0.5) + // upper pixel = (int) (p + 0.5) + + // XXX use 'floor' function so -0.2 becomes -1 + int pm = floor (p - 0.5); + int pn = floor (p + 0.5); + + // if pm = -1, extrapolate from p[0] and p[1]: + // if pn = Npix, extrapolate from p[-1] and p[-2]: + if (pm == -1) { pm ++; pn ++; } + if (pn == Npix) { pm --; pn --; } + + // 5.6 - 5 - 0.5 = 0.1 + // 0.4 - 0 - 0.5 = -0.1 + float f = p - pm - 0.5; + float Pval = source->profile->data.F32[pm]*(1.0 - f) + source->profile->data.F32[pn]*f; + + float value = Pval * exp (-0.5 * PS_SQR(q) / PS_SQR(sigma)); + image->data.F32[iy][ix] += value; + } + } + + // dy = p[1] - p[0] + // dx = 1 + // y' = p[0] + f*dy = p[0] + f*p[1] - f*p[0] + + // x = 0.3 + // x - 0.5 = -0.2 + // pm = -1 + // pm -> 0, pn -> 1 + // f = 0.3 + 1 - 0.5 = 0.8 + + // dy = p[1] - p[0] + // dx = 1 + // y = p[0] + + return true; +} + +// extract source from image given centroid and posangle +// nearly there, but off by 1 in the x-direction +psImage *pmSourceTraceGet (psImage *image, pmSourceTrace *source) { + + if (!image) return NULL; + if (!source) return NULL; + + // given a trace with a specified centroid and posangle, determine sigma and the profile + + float CS = dCOS(source->posAngle); + float SN = dSIN(source->posAngle); + + // int Nx = image->numCols; + // int Ny = image->numRows; + + // first, extract the pixels from the image to the trace image in P,Q pixel frame + psImage *window = psImageAlloc(source->dP, source->dQ, PS_TYPE_F32); + psImageInit (window, 0.0); + + // XXX note these should be the same as source->dP,dQ + int Np = window->numCols; + int Nq = window->numRows; + + // centroid corresponds to center of pixel 0, + for (int iq = 0; iq < Nq; iq++) { + float dq = floor(iq - 0.5*Nq) + 0.5; + for (int ip = 0; ip < Np; ip++) { + float dp = ip + 0.5; + + // transform (dp,dq) to (x,y) where p is along the profile an q is cross profile + // (dp,dq) = (0,0) corresponds to src->(x,y) + // center line of the trace is midpoint of the window + + float dx = CS*dp - SN*dq; + float dy = CS*dq + SN*dp; + + // pixel coordinate in image coords + float fx = dx + source->x; + float fy = dy + source->y; + + int ixm = floor(fx); + int ixn = ixm + 1; + + int iym = floor(fy); + int iyn = iym + 1; + + float dfx = fx - ixm; + float dfy = fy - iym; + + float V00 = image->data.F32[iym][ixm]; + float V01 = image->data.F32[iym][ixn]; + float V10 = image->data.F32[iyn][ixm]; + float V11 = image->data.F32[iyn][ixn]; + + float V0 = (1.0 - dfx)*V00 + dfx*V01; + float V1 = (1.0 - dfx)*V10 + dfx*V11; + + window->data.F32[iq][ip] = (1.0 - dfy)*V0 + dfy*V1; + } + } + return window; +} + +// XXX defined in pmSourceTraceFitGauss.c +bool myGaussFit (psVector *xValues, psVector *yValues); + +// given window with extracted source trace, calculate the profile and sigma +psVector *pmSourceTraceExtractProfile (psImage *window, pmSourceTrace *source) { + + if (!window) return false; + if (!source) return false; + + // XXX note these should be the same as source->dP,dQ + int Np = window->numCols; // x-direction is along the profile + int Nq = window->numRows; // y-direction is cross-profile + + // is the window background-subtracted? must be, otherwise I need more area + + // Each column is a measurement of the cross-profile PSF, but each has a different + // normalization. Measure the integrated flux across the profile and divide by that + // as a measurement of the normalization. Adjust the weight of each pixel based on the + // normalization. Sum across columns for each pixels using + + // sums contains the normalization for each column + psVector *sums = psVectorAlloc (Np, PS_TYPE_F32); + psVectorInit (sums, 0.0); + + for (int ip = 0; ip < Np; ip++) { + for (int iq = 0; iq < Nq; iq++) { + sums->data.F32[ip] += window->data.F32[iq][ip]; + } + } + + // crosscut will contain the cross-profile PSF + psVector *crosscut = psVectorAlloc (Nq, PS_TYPE_F32); + psVectorInit (crosscut, 0.0); + + psVector *crosswgt = psVectorAlloc (Nq, PS_TYPE_F32); + psVectorInit (crosswgt, 0.0); + + for (int iq = 0; iq < Nq; iq++) { + for (int ip = 0; ip < Np; ip++) { + // this gives each column effectively equal weight, which seems wrong + // crosscut->data.F32[iq] += window->data.F32[iq][ip] / sums->data.F32[ip]; + // crosswgt->data.F32[iq] += 1.0; + + // this weights each column by the S/N of that column + if (fabs(sums->data.F32[ip]) < 1) continue; + crosscut->data.F32[iq] += window->data.F32[iq][ip]; + crosswgt->data.F32[iq] += 1.0 / sums->data.F32[ip]; + } + } + + for (int iq = 0; iq < Nq; iq++) { + crosscut->data.F32[iq] /= crosswgt->data.F32[iq]; + } + + psVector *xCoord = psVectorAlloc(crosscut->n, PS_TYPE_F32); + for (int i = 0; i < xCoord->n; i++) { xCoord->data.F32[i] = i; } + + myGaussFit (xCoord, crosscut); + + psFree (crosswgt); + psFree (xCoord); + return crosscut; +} + Index: /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSourceTrace.h =================================================================== --- /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSourceTrace.h (revision 42955) +++ /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSourceTrace.h (revision 42955) @@ -0,0 +1,39 @@ +/* @file pmSourceTrace.h + * + * @author EAM, IfA + * + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2014-03-20 02:31:25 $ + * Copyright 2014 IfA, University of Hawaii + */ + +# ifndef PM_SOURCE_TRACE_H +# define PM_SOURCE_TRACE_H + +// structure to describe a model for a source which is a linear trail +// convolved in the cross direction by the PSF. In one direction, the +// source flux is described by a profile with flux vs offset. The entire +// profile has a specific position angle +typedef struct { + psVector *profile; // normalized profile? + + /// XXX for the initial implementation, assume a Gaussian cross profile + psF32 Io; // is this needed (integrate into the profile?) + psF32 sigma; // PSF Gaussian sigma + psF32 posAngle; // position angle of the trace (CCW from x-axis) + psF32 x; // reference pixel in image coordinates + psF32 y; // reference pixel in image coordinates + + // should the profile define the full window in that direction? + // trace can be assumed centered on the profile + + int dP; // window size in profile direction + int dQ; // window size in cross-profile direction +} pmSourceTrace; + +pmSourceTrace *pmSourceTraceAlloc(psF32 x, psF32 y, psF32 posAngle, psF32 sigma); +bool pmSourceTraceAdd (psImage *image, pmSourceTrace *source); +psImage *pmSourceTraceGet (psImage *image, pmSourceTrace *source); +psVector *pmSourceTraceExtractProfile (psImage *window, pmSourceTrace *source); + +# endif /* PM_SOURCE_TRACE_H */ Index: /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSourceTraceFitGauss.c =================================================================== --- /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSourceTraceFitGauss.c (revision 42955) +++ /branches/eam_branches/ipp-20250626/psModules/src/objects/pmSourceTraceFitGauss.c (revision 42955) @@ -0,0 +1,141 @@ +/** @file pmSourceTrace.c + * @author EAM, IfA + * @date $Date: 2008-10-08 21:53:08 $ + * Copyright 2025 IfA, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +// #include "pmSourceTrace.h" + +#define NUM_ITERATIONS 100 +#define ERR_TOL 1e-6 + +float myFindThreshold (psVector *x, psVector *y, float value); + +// mean = params->data.F32[0] +// sigma = params->data.F32[1] +// norm = params->data.F32[2] + +// y = p2 * exp (-0.5*(x - p0)^2 / p1^2) +// speed this up by making it explicitly 1D? +// (would need to create new psMinimize func) +float myGauss(const psVector *params, psVector *x) +{ + // XXX speed up by fitting for sigma^2 ? + float z = (x->data.F32[0] - params->data.F32[0]) / params->data.F32[1]; + float r = expf (-0.5 * PS_SQR (z)); + float value = params->data.F32[2] * r; + return value; +} + +void myGaussDeriv(psVector *deriv, const psVector *params, psVector *x) +{ + float z = (x->data.F32[0] - params->data.F32[0]) / params->data.F32[1]; + float r = expf (-0.5 * PS_SQR (z)); + + deriv->data.F32[0] = params->data.F32[2]*r*z/params->data.F32[1]; + deriv->data.F32[1] = params->data.F32[2]*r*z*z/params->data.F32[1]; + deriv->data.F32[2] = r; +} + +psF32 myGaussFitFunc(psVector *deriv, + psVector *params, + psVector *x) +{ + if ((deriv == NULL) || (params == NULL) || (x == NULL)) { + psError(PS_ERR_UNKNOWN, true, "deriv or params or x is NULL.\n"); + } + + myGaussDeriv(deriv, params, x); + return myGauss(params, x); +} + +bool myGaussFit (psVector *xValues, psVector *yValues) { + + psMinimization *min = psMinimizationAlloc(NUM_ITERATIONS, ERR_TOL, 100); + + int Npts = xValues->n; + // XXX check yValues->n + + psArray *coords = psArrayAlloc(Npts); + psVector *values = psVectorAlloc(Npts, PS_TYPE_F32); + psVector *errors = psVectorAlloc(Npts, PS_TYPE_F32); + psVector *params = psVectorAlloc(3, PS_TYPE_F32); + + // XXX set the guess params + psVector *cumulate = psVectorAlloc(Npts, PS_TYPE_F32); + psVectorInit (cumulate, 0.0); + + float max = yValues->data.F32[0]; + float sum = 0.0; + for (int i = 0; i < Npts; i++) { + sum += yValues->data.F32[i]; + cumulate->data.F32[i] = sum; + max = PS_MAX(max, yValues->data.F32[i]); + } + for (int i = 0; i < Npts; i++) { + cumulate->data.F32[i] /= sum; + } + + // where does cumulate exceed 0.16 and 0.84 + float f16 = myFindThreshold (xValues, cumulate, 0.16); + float f50 = myFindThreshold (xValues, cumulate, 0.50); + float f84 = myFindThreshold (xValues, cumulate, 0.84); + + fprintf (stderr, "f16, f50, f84: %f, %f, %f\n", f16, f50, f84); + + params->data.F32[0] = f50; + params->data.F32[1] = (f84 - f16) / 2.0; + params->data.F32[2] = max; + + fprintf (stderr, "p0, p1, p2: %f, %f, %f\n", params->data.F32[0], params->data.F32[1], params->data.F32[2]); + + for (long i = 0; i < Npts; i++) { + psVector *x = psVectorAlloc(1, PS_TYPE_F32); + x->data.F32[0] = xValues->data.F32[i]; + coords->data[i] = x; + + values->data.F32[i] = yValues->data.F32[i]; + errors->data.F32[i] = 0.1; // use const or poisson? + } + + bool myResult = psMinimizeLMChi2(min, NULL, params, NULL, coords, values, errors, (psMinimizeLMChi2Func)myGaussFitFunc); + printf("Minimisation took %d iterations\n", min->iter); + printf("chi^2 at the minimum is %.3g\n", min->value); + for (long i = 0; i < params->n; i++) { + fprintf (stderr, "Parameter %ld at the minimum is %.3f\n", i, params->data.F32[i]); + } + + psFree(min); + psFree(params); + psFree(coords); + psFree(values); + psFree(errors); + return myResult; +} + +float myFindThreshold (psVector *x, psVector *y, float value) { + + int i; + for (i = 0; (i < x->n) && (y->data.F32[i] < value); i++); + + // linear interpolation between i and i-1 + float dy = value - y->data.F32[i-1]; + + float dxdy = (x->data.F32[i] - x->data.F32[i-1]) / (y->data.F32[i] - y->data.F32[i-1]); + + float xValue = x->data.F32[i-1] + dxdy * dy; + + return xValue; + +} + +