Index: trunk/Ohana/src/opihi/mana/Makefile =================================================================== --- trunk/Ohana/src/opihi/mana/Makefile (revision 41158) +++ trunk/Ohana/src/opihi/mana/Makefile (revision 41159) @@ -26,8 +26,8 @@ $(SRC)/deimos_getobj.$(ARCH).o \ $(SRC)/deimos_fitobj.$(ARCH).o \ -$(SRC)/deimos_fitalt.$(ARCH).o \ -$(SRC)/deimos_mkalt.$(ARCH).o \ $(SRC)/deimos_mkmodel.$(ARCH).o \ $(SRC)/deimos_arclines.$(ARCH).o \ +$(SRC)/deimos_fitarc.$(ARCH).o \ +$(SRC)/deimos_fitprofile.$(ARCH).o \ $(SRC)/findrowpeaks.$(ARCH).o Index: trunk/Ohana/src/opihi/mana/deimos.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos.c (revision 41158) +++ trunk/Ohana/src/opihi/mana/deimos.c (revision 41159) @@ -4,19 +4,23 @@ int deimos_fitslit (int argc, char **argv); int deimos_mkobj (int argc, char **argv); -int deimos_mkalt (int argc, char **argv); int deimos_getobj (int argc, char **argv); +// int deimos_getalt (int argc, char **argv); int deimos_fitobj (int argc, char **argv); -int deimos_fitalt (int argc, char **argv); +// int deimos_fitalt (int argc, char **argv); int deimos_arclines (int argc, char **argv); +int deimos_fitarc (int argc, char **argv); +int deimos_fitprofile (int argc, char **argv); static Command deimos_commands[] = { - {1, "fitobj", deimos_fitobj, "fit for object parameters"}, - {1, "fitalt", deimos_fitalt, "fit for object parameters using LMM"}, - {1, "getobj", deimos_getobj, "determine crude object parameters"}, - {1, "mkobj", deimos_mkobj, "make a full object image"}, - {1, "mkalt", deimos_mkalt, "make a full object image (uses deimos_make_object, for a test)"}, - {1, "mkslit", deimos_mkslit, "make a slit image"}, - {1, "fitslit", deimos_fitslit, "fit slit image to observed slit flux"}, - {1, "arclines", deimos_arclines, "detect arclines using LSF and STILT"}, + {1, "fitobj", deimos_fitobj, "fit for object parameters using LMM"}, + {1, "fitalt", deimos_fitobj, "fit for object parameters using LMM"}, + {1, "getobj", deimos_getobj, "determine crude object parameters"}, + {1, "getalt", deimos_getobj, "determine crude object parameters"}, + {1, "mkobj", deimos_mkobj, "make a full object image"}, + {1, "mkslit", deimos_mkslit, "make a slit image"}, + {1, "fitslit", deimos_fitslit, "fit slit image to observed slit flux"}, + {1, "arclines", deimos_arclines, "detect arclines using LSF and STILT"}, + {1, "fitarc", deimos_fitarc, "fit arclamp lines"}, + {1, "fitprofile", deimos_fitprofile, "fit slit profile"}, }; Index: trunk/Ohana/src/opihi/mana/deimos_arclines.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_arclines.c (revision 41158) +++ trunk/Ohana/src/opihi/mana/deimos_arclines.c (revision 41159) @@ -53,7 +53,14 @@ deimos_make_LSF_kernel (LSF, stilt, Nx); + // FILE *Fout = NULL; + // loop over the rows for (int iy = 0; iy < Ny; iy++) { float Fsum = 0; + float Ksum = 0; + + // if (iy == 5802) Fout = fopen ("test.5802.dat", "w"); + // if (iy == 5815) Fout = fopen ("test.5815.dat", "w"); + // cross-correlation of buffer values and kernel values for (int ix = 0; ix < Nx; ix++) { @@ -62,7 +69,16 @@ if ((iy + ky) < 0) continue; if ((iy + ky) >= Ny) continue; + if (!isfinite(Fin[ix + (iy + ky)*Nx])) continue; + if (!isfinite(kernel[ix + ko*Nx])) continue; Fsum += Fin[ix + (iy + ky)*Nx] * kernel[ix + ko*Nx]; + Ksum += kernel[ix + ko*Nx]; + // if ((iy == 5815) || (iy == 5802)) { + // fprintf (Fout, "%d %d %d : %f %f\n", iy, ix, ky, Fin[ix + (iy + ky)*Nx], kernel[ix + ko*Nx]); + // } } } + + // if ((iy == 5815) || (iy == 5802)) fclose (Fout); + coord->elements.Flt[iy] = iy; flux->elements.Flt[iy] = Fsum; Index: trunk/Ohana/src/opihi/mana/deimos_fitalt.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_fitalt.c (revision 41158) +++ (revision ) @@ -1,447 +1,0 @@ -# include "data.h" -# include "deimos.h" - -# if (0) -int deimos_fitalt (int argc, char **argv) { - OHANA_UNUSED_PARAM(argc); - OHANA_UNUSED_PARAM(argv); - return FALSE; -} - -# else - -/* - this is starting to work OK. some improvements to make - * use Gaussdev to sample - * do not scale down range (or user-set scale-down) - */ - -// internal functions to fitobj -static float deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts); -double deimos_LMM_update (float *data, float *model_ref, float *model_obj, float *model_sky, float *model_bck, float *dObj, float *dSky, float *dBck, int Nx, int Ny, int row, float lambda); -static opihi_flt *deimos_make_test (opihi_flt *guess, float sigma, int row, int Nrow, float *dValue); -opihi_flt *deimos_subvector (opihi_flt *fullVec, int row, int Nrow); - -int deimos_fitalt (int argc, char **argv) { - - // input parameters: - // * buffer : 2D image of slit (full frame or cutout?) - // * trace : spline fit of slit central x pos vs y-coord - // * profile : slit window profile (vector) - // * PSF : point-spread function vector (flux normalized, x-dir) - // * stilt : slit tilt response : 2D kernel? - - // in-out parameters: - // * obj : vector of object flux vs y-coord (starting guess and result) - // * sky : vector of local sky signal vs y-coord - // * bck : vector of extra-slit background flux vs y-coord - - int N; - - Spline *trace = NULL; - Vector *profile = NULL; - - Vector *obj = NULL; - Vector *sky = NULL; - Vector *bck = NULL; - - Vector *psf = NULL; - - Buffer *buffer = NULL; - - float stilt = 0.0; // angle of the slit - ohana_gaussdev_init (); - - int DO_PLUS = FALSE; - if ((N = get_argument (argc, argv, "-plus"))) { - remove_argument (N, &argc, argv); - DO_PLUS = TRUE; - } - - int VERBOSE = FALSE; - if ((N = get_argument (argc, argv, "-v"))) { - remove_argument (N, &argc, argv); - VERBOSE = TRUE; - } - - // a crude noise model valid for both the guess vectors and the - // data buffer: - float gain = 1.0; - if ((N = get_argument (argc, argv, "-gain"))) { - remove_argument (N, &argc, argv); - gain = atof(argv[N]); - remove_argument (N, &argc, argv); - } - float noise = 0.0; - if ((N = get_argument (argc, argv, "-noise"))) { - remove_argument (N, &argc, argv); - noise = atof(argv[N]); - remove_argument (N, &argc, argv); - } - float nsigma = 2.0; - if ((N = get_argument (argc, argv, "-nsigma"))) { - remove_argument (N, &argc, argv); - nsigma = atof(argv[N]); - remove_argument (N, &argc, argv); - } - float lambda = 1.0; - if ((N = get_argument (argc, argv, "-lambda"))) { - remove_argument (N, &argc, argv); - lambda = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - int Niter = 25; - if ((N = get_argument (argc, argv, "-iter"))) { - remove_argument (N, &argc, argv); - Niter = atoi(argv[N]); - remove_argument (N, &argc, argv); - } - int Nrow = 5; - if ((N = get_argument (argc, argv, "-row"))) { - remove_argument (N, &argc, argv); - Nrow = atoi(argv[N]); - remove_argument (N, &argc, argv); - } - - // Input parameters: - if ((N = get_argument (argc, argv, "-trace"))) { - remove_argument (N, &argc, argv); - if ((trace = FindSpline (argv[N])) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } else { goto usage; } - if ((N = get_argument (argc, argv, "-profile"))) { - remove_argument (N, &argc, argv); - if ((profile = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } else { goto usage; } - if ((N = get_argument (argc, argv, "-psf"))) { - remove_argument (N, &argc, argv); - if ((psf = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } else { goto usage; } - if ((N = get_argument (argc, argv, "-stilt"))) { - remove_argument (N, &argc, argv); - stilt = atof (argv[N]); - remove_argument (N, &argc, argv); - } else { goto usage; } - - // In-Out parameters: - if ((N = get_argument (argc, argv, "-object"))) { - remove_argument (N, &argc, argv); - if ((obj = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } else { goto usage; } - if ((N = get_argument (argc, argv, "-sky"))) { - remove_argument (N, &argc, argv); - if ((sky = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } else { goto usage; } - if ((N = get_argument (argc, argv, "-backgnd"))) { - remove_argument (N, &argc, argv); - if ((bck = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); - remove_argument (N, &argc, argv); - } else { goto usage; } - - if (argc != 3) goto usage; - - // buffer is the full image, Xref is reference coordinate of the profile in the buffer - if ((buffer = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); - - // profile and PSF are defined with central reference pixel mapped to Xref (+ trace) - if (profile->Nelements % 2 == 0) { - gprint (GP_ERR, "slit profile vector must have an odd number of pixels\n"); - return FALSE; - } - if (psf->Nelements % 2 == 0) { - gprint (GP_ERR, "PSF vector must have an odd number of pixels\n"); - return FALSE; - } - - // observed buffer size - int Nx = buffer[0].matrix.Naxis[0]; - int Ny = buffer[0].matrix.Naxis[1]; - float *bufVal = (float *) buffer->matrix.buffer; - - // obj, sky, bck must be consistent with data - if (Ny != obj->Nelements) { - gprint (GP_ERR, "inconsistent wavelength scales (object)\n"); - return FALSE; - } - if (Ny != sky->Nelements) { - gprint (GP_ERR, "inconsistent wavelength scales (sky)\n"); - return FALSE; - } - if (Ny != bck->Nelements) { - gprint (GP_ERR, "inconsistent wavelength scales (backgnd)\n"); - return FALSE; - } - - // for the cross-dispersion reference pixel, use the user value or set to Nx/2 if < 0 - deimos_set_cross_ref (atoi(argv[2]), Nx); - deimos_make_kernel (stilt, Nx); - - Vector *objMin = NULL; - Vector *skyMin = NULL; - Vector *bckMin = NULL; - - if ((objMin = SelectVector ("objMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (objMin, OPIHI_FLT, Ny); - if ((skyMin = SelectVector ("skyMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (skyMin, OPIHI_FLT, Ny); - if ((bckMin = SelectVector ("bckMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bckMin, OPIHI_FLT, Ny); - - // the functions below accept an opihi_flt array for object, sky, background - opihi_flt *objVal = objMin->elements.Flt; - opihi_flt *skyVal = skyMin->elements.Flt; - opihi_flt *bckVal = bckMin->elements.Flt; - - // we make a copy of the input guess vectors and minimize them - for (int iy = 0; iy < Ny; iy++) { - objVal[iy] = obj->elements.Flt[iy]; - skyVal[iy] = sky->elements.Flt[iy]; - bckVal[iy] = bck->elements.Flt[iy]; - } - - float noiseVar = SQ(noise); - ALLOCATE_PTR (objNoise, float, Ny); - ALLOCATE_PTR (skyNoise, float, Ny); - ALLOCATE_PTR (bckNoise, float, Ny); - - // generate noise vectors - for (int iy = 0; iy < Ny; iy++) { - objNoise[iy] = nsigma*sqrt(hypot(fabs(objVal[iy]) / gain, noiseVar)); - skyNoise[iy] = nsigma*sqrt(hypot(fabs(skyVal[iy]) / gain, noiseVar)); - bckNoise[iy] = nsigma*sqrt(hypot(fabs(bckVal[iy]) / gain, noiseVar)); - } - - // obj, sky, bck are initial guesses. Also use these values to define a range for the initial guesses. - - // XXX worry about last segment (< Nrow rows) - for (int iter = 0; iter < Niter; iter++) { - - int dRow = Nrow/2; - - // use scaled-noise to define test ranges - // we are going to try each row using a window Nrow wide to measure chisq: - for (int row = dRow; row < Ny - dRow - 1; row++) { - - // I am using a numerical equivalent of Levenberg-Marquardt, but only for a single - // p_m[i] = (obj,sky,bck) triplet [m is the index of obj,sky,bck; i is the row - // index] at a time. I start with a current best guess set of (obj,sky,bck). I - // need to calculate the model at this location, and then at offset locations p_m + - // dp_m. Note that I do not need to calculate chisq for each of these to generate - // the elements of the LMM equations: - - int row_0 = row - dRow; // start row of the subset test range - - // Nrow = 5, dRow = 2, Ny = 100, last start is Ny - 2 - 1 = 97 - // last row_0 is 95 - - // generate subset vectors for the range **centered** on row - opihi_flt *objCurr = deimos_subvector (objVal, row_0, Nrow); - opihi_flt *skyCurr = deimos_subvector (skyVal, row_0, Nrow); - opihi_flt *bckCurr = deimos_subvector (bckVal, row_0, Nrow); - - // current vector value: - // for deimos_make_model, row is the starting point of the subimage - float *model_ref = deimos_make_model (objCurr, skyCurr, bckCurr, psf, profile, trace, Nx, Nrow, row_0); - - // in the test functions below, dRow is the central pixel of the subset, e.g., obj[row] => objCurr[dRow] - - float dObj = 0.0, dSky = 0.0, dBck = 0.0; - - // delta obj (save dObj) - opihi_flt *objtest = deimos_make_test (objCurr, objNoise[row], dRow, Nrow, &dObj); - float *model_obj = deimos_make_model (objtest, skyCurr, bckCurr, psf, profile, trace, Nx, Nrow, row_0); - - // delta sky (save dSky) - opihi_flt *skytest = deimos_make_test (skyCurr, skyNoise[row], dRow, Nrow, &dSky); - float *model_sky = deimos_make_model (objCurr, skytest, bckCurr, psf, profile, trace, Nx, Nrow, row_0); - - // delta bck (save dBck) - opihi_flt *bcktest = deimos_make_test (bckCurr, bckNoise[row], dRow, Nrow, &dBck); - float *model_bck = deimos_make_model (objCurr, skyCurr, bcktest, psf, profile, trace, Nx, Nrow, row_0); - - // dObj, dSky, dBck carry the delta used derivatives in and the delta applied to the values out - float chisq = deimos_LMM_update (bufVal, model_ref, model_obj, model_sky, model_bck, &dObj, &dSky, &dBck, Nx, Nrow, row_0, lambda); - - if (1) { - if (DO_PLUS) { - objVal[row] += dObj; - skyVal[row] += dSky; - bckVal[row] += dBck; - } else { - objVal[row] -= dObj; - skyVal[row] -= dSky; - bckVal[row] -= dBck; - } - } - - free (objCurr); - free (skyCurr); - free (bckCurr); - - free (model_ref); - - free (model_obj); - free (model_sky); - free (model_bck); - - free (objtest); - free (skytest); - free (bcktest); - - if (VERBOSE) fprintf (stderr, "chisq: %f (%f, %f, %f)\n", chisq, dObj, dSky, dBck); - } - - float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, trace, Nx, Ny, 0); - - int Npts = 0; - float chisq_full = deimos_get_chisq (bufVal, model_full, Nx, Ny, 0, &Npts); - - fprintf (stderr, "** Full chisq: %f (%d) : %f)\n", chisq_full, Npts, chisq_full / (1.0 * Npts)); - } - - return TRUE; - - usage: - gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n"); - return FALSE; -} - -/****************** fitobj Support Functions *******************/ - -// we generate a new vector with a single element (row) modified based on sigma -opihi_flt *deimos_subvector (opihi_flt *fullVec, int row, int Nrow) { - - ALLOCATE_PTR (value, opihi_flt, Nrow); - - for (int i = 0; i < Nrow; i++) { - value[i] = fullVec[i + row]; - } - - return value; -} - - -// we generate a new vector with a single element (row) modified based on sigma -static opihi_flt *deimos_make_test (opihi_flt *guess, float sigma, int row, int Nrow, float *dValue) { - - // if fullInput is TRUE, the guess vector runs from 0 to Ny (full wavelength range) - // if fullInput is FALSE, the guess vector runs from row to row + Nrow (subset range) - // sigma is valid for the element at row - - ALLOCATE_PTR (value, opihi_flt, Nrow); - - for (int i = 0; i < Nrow; i++) { - value[i] = guess[i]; - } - - value[row] = guess[row] + 0.1*sigma; - - *dValue = 0.1*sigma; - return value; -} - - -double deimos_LMM_update (float *data, float *model_ref, - float *model_obj, float *model_sky, float *model_bck, - float *dObj, float *dSky, float *dBck, int Nx, int Ny, int row, float lambda) { - - // we are going to loop over the images calculating the following value: - - // W = 1 / sigma^2, leave as 1.0 for now - // chisq = sum W * (data - model_ref)^2 - // dF_X = sum W * (data - model_ref) * (model_ref - model_X) / dX - // d2F_XY = sum W * (model_ref - model_X) * (model_ref - model_Y) / (dX * dY) - - double chisq = 0.0; - double dF_obj = 0.0, dF_sky = 0.0, dF_bck = 0.0; - double d2F_obj_obj = 0.0, d2F_obj_sky = 0.0, d2F_obj_bck = 0.0, d2F_sky_sky = 0.0, d2F_sky_bck = 0.0, d2F_bck_bck = 0.0; - - float W = 1.0; - - for (int iy = 0; iy < Ny; iy++) { - for (int ix = 0; ix < Nx; ix++) { - int pix = ix + iy*Nx; - int pix_data = ix + (iy + row)*Nx; - - // if any of the pixel values are NAN, skip the point: - if (!isfinite(data[pix_data])) continue; - if (!isfinite(model_ref[pix])) continue; - if (!isfinite(model_obj[pix])) continue; - if (!isfinite(model_sky[pix])) continue; - if (!isfinite(model_bck[pix])) continue; - - float df = W * (data[pix_data] - model_ref[pix]); - - chisq += SQ(df); - - float df_obj = (model_ref[pix] - model_obj[pix]) / *dObj; - float df_sky = (model_ref[pix] - model_sky[pix]) / *dSky; - float df_bck = (model_ref[pix] - model_bck[pix]) / *dBck; - - dF_obj += df*df_obj; - dF_sky += df*df_sky; - dF_bck += df*df_bck; - - d2F_obj_obj += W * SQ(df_obj); - d2F_sky_sky += W * SQ(df_sky); - d2F_bck_bck += W * SQ(df_bck); - - d2F_obj_sky += W * df_obj * df_sky; - d2F_obj_bck += W * df_obj * df_bck; - d2F_sky_bck += W * df_sky * df_bck; - } - } - - // these are the elements of Ax = B - ALLOCATE_PTR (A, double *, 3); - ALLOCATE_PTR (B, double *, 3); - for (int i = 0; i < 3; i++) { - ALLOCATE (A[i], double, 3); - ALLOCATE (B[i], double, 1); - } - - A[0][0] = d2F_obj_obj * (1.0 + lambda); - A[1][1] = d2F_sky_sky * (1.0 + lambda); - A[2][2] = d2F_bck_bck * (1.0 + lambda); - - A[0][1] = A[1][0] = d2F_obj_sky; - A[1][2] = A[2][1] = d2F_sky_bck; - A[0][2] = A[2][0] = d2F_obj_bck; - - B[0][0] = dF_obj; - B[1][0] = dF_sky; - B[2][0] = dF_bck; - - dgaussjordan (A, B, 3, 1); - - *dObj = B[0][0]; - *dSky = B[1][0]; - *dBck = B[2][0]; - - return chisq; -} - -static float deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts) { - - int npts = 0; - float chisq = 0; - - for (int iy = 0; iy < Ny; iy++) { - - for (int ix = 0; ix < Nx; ix++) { - - int pix_buffer = ix + (iy + row)*Nx; - int pix_model = ix + iy*Nx; - - if (!isfinite(buffer[pix_buffer])) continue; - if (!isfinite(model[pix_model])) continue; - chisq += SQ(buffer[pix_buffer] - model[pix_model]); - npts ++; - } - } - *Npts = npts; - return chisq; -} -# endif Index: trunk/Ohana/src/opihi/mana/deimos_fitarc.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_fitarc.c (revision 41159) +++ trunk/Ohana/src/opihi/mana/deimos_fitarc.c (revision 41159) @@ -0,0 +1,237 @@ +# include "data.h" +# define NPAR 4 + +typedef opihi_flt FitFunc (opihi_flt, opihi_flt *, int, opihi_flt *); + +opihi_flt fgaussOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar); + +void gaussian_fit_irls (opihi_flt *xprofile, opihi_flt *fprofile, opihi_flt *wprofile, opihi_flt *oprofile, int Npts, opihi_flt *par, int Npar); +void gaussian_fit (opihi_flt *xprofile, opihi_flt *fprofile, opihi_flt *wprofile, opihi_flt *oprofile, int Npts, opihi_flt *par, int Npar, FitFunc *FUNC); + +double deimos_weight_cauchy (double x); + +// use these globals to carry into sigmoid_fit +static int VERBOSE = FALSE; +static int QUIET = FALSE; +static int MaxIterations = 10; + +/* fitarc buffer Xo Yo Nx Ny + + * for each col, fit the y-dir profile with a gauss + * fit the C0 values to a line to get the angle + */ + +int deimos_fitarc (int argc, char **argv) { + + // input parameters: + // * xprofile : extracted slit window coordinates (vector) + // * xprofile : extracted slit window profile (vector) + // * Xs, Xe : scale coordinates of window boundaries + + // output parameters: + // * sigmoid_left_center, etc : coordinates of window boundaries + + int N; + Buffer *buffer = NULL; + + QUIET = FALSE; + if ((N = get_argument (argc, argv, "-q"))) { + QUIET = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-quiet"))) { + QUIET = TRUE; + remove_argument (N, &argc, argv); + } + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + MaxIterations = 10; + if ((N = get_argument (argc, argv, "-max-iterations"))) { + remove_argument (N, &argc, argv); + MaxIterations = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + Buffer *residbuf = NULL; + if ((N = get_argument (argc, argv, "-resid"))) { + remove_argument (N, &argc, argv); + residbuf = SelectBuffer (argv[N], ANYBUFFER, TRUE); + if (!residbuf) goto usage; + remove_argument (N, &argc, argv); + } + + Buffer *varbuf = NULL; + if ((N = get_argument (argc, argv, "-varbuf"))) { + remove_argument (N, &argc, argv); + varbuf = SelectBuffer (argv[N], OLDBUFFER, TRUE); + if (!varbuf) goto usage; + remove_argument (N, &argc, argv); + } + + if (argc != 6) goto usage; + if ((buffer = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + + int nx = buffer[0].matrix.Naxis[0]; + int ny = buffer[0].matrix.Naxis[1]; + + // supplied guess for boundaries + int Xs = atof(argv[2]); + int Ys = atof(argv[3]); + int Nx = atof(argv[4]); + int Ny = atof(argv[5]); + + if (residbuf) { + ResetBuffer (residbuf, nx, ny, -32, 0.0, 1.0); + } + + // validate Xs, Ys, Nx, Ny + // Xs + Nx <= nx + + // loop over the cols + + ALLOCATE_PTR (xcol, double, Ny); + ALLOCATE_PTR (ycol, double, Ny); + ALLOCATE_PTR (ocol, double, Ny); // fitted value + + float *bufVal = (float *) buffer->matrix.buffer; + + for (int ix = Xs; ix < Xs + Nx; ix++) { + + // loop over the pixels in this column to get the x and y vectors + // find the peak at the same time + + double xmax = -HUGE_VAL; + double ymax = -HUGE_VAL; + int Nc = 0; + + for (int iy = 0; iy < Ny; iy++) { + + float yval = bufVal[ix + (iy + Ys)*Nx]; + if (!isfinite(yval)) continue; + + xcol[Nc] = iy + Ys; + ycol[Nc] = yval; + + if (yval > ymax) { + ymax = yval; + xmax = iy + Ys; + } + Nc ++; + } + + // now do a Gauss fit given C0 = xmax, C1 = 2, C2 = ymax, C3 = 0 + + // temporary storage vectors + opihi_flt par[NPAR]; + par[0] = xmax; + par[1] = 2; + par[2] = ymax; + par[3] = 0; + gaussian_fit_irls (xcol, ycol, NULL, ocol, Nc, par, NPAR); + + fprintf (stderr, "%d : %f %f : %f %f\n", ix, par[0], par[2], par[1], par[3]); + set_variable ("C0", par[0]); + set_variable ("C1", par[1]); + set_variable ("C2", par[2]); + set_variable ("C3", par[3]); + + if (residbuf) { + float *residVal = (float *) residbuf->matrix.buffer; + for (int iy = 0; iy < Ny; iy++) { + float yval = fgaussOD (iy + Ys, par, NPAR, NULL); + // residVal[ix + (iy + Ys)*Nx] = bufVal[ix + (iy + Ys)*Nx] - yval; + residVal[ix + (iy + Ys)*Nx] = yval; + } + } + } + + return TRUE; + + usage: + gprint (GP_ERR, "USAGE: deimos fitarc buff Xo Yo Nx Ny\n"); + gprint (GP_ERR, " [-resid buffer] : output residual image\n"); + gprint (GP_ERR, " [-varbuf buffer] : input variance image\n"); + return FALSE; +} + +# define FIT_TOLERANCE 1e-4 +# define FLT_TOLERANCE 1e-6 + +// xprofile : x-coordinate +// fprofile : flux (counts) +// wprofile : weight (true poisson error) +// oprofile : output profile (not used?) +void gaussian_fit_irls (opihi_flt *xprofile, opihi_flt *fprofile, opihi_flt *wprofile, opihi_flt *oprofile, int Npts, opihi_flt *par, int Npar) { + + // buffer to store result from each iteration + ALLOCATE_PTR (par_old, opihi_flt, Npar); + + // define initial weights as 1.0 + ALLOCATE_PTR (wt, opihi_flt, Npts); + ALLOCATE_PTR (wtemp, opihi_flt, Npts); + for (int i = 0; i < Npts; i++) { + wtemp[i] = wprofile ? (fabs(wprofile[i]) > 1e-5 ? wprofile[i] : 1e-5) : 1.0; + wt[i] = 1.0 / SQ(wtemp[i]); + } + + // initial fit with ordinary least-squares + // otemp is the fitted value of the function at xtemp + gaussian_fit (xprofile, fprofile, wt, NULL, Npts, par, Npar, fgaussOD); + + int converged = FALSE; + for (int iterations = 0; !converged && (iterations < MaxIterations); iterations++) { + + for (int i = 0; i < Npts; i++) { + wt[i] = deimos_weight_cauchy ((fprofile[i] - oprofile[i]) / wtemp[i]) / SQ(wtemp[i]); + } + + gaussian_fit (xprofile, fprofile, wt, oprofile, Npts, par, Npar, fgaussOD); + + // save this solution + for (int i = 0; i < Npar; i++) { + par_old[i] = par[i]; + } + + converged = TRUE; + for (int i = 0; i < Npts; i++) { + if ((fabs(par[i] - par_old[i]) > FIT_TOLERANCE * fabs(par[i])) && + (fabs(par[i] - par_old[i]) > FLT_TOLERANCE)) + converged = FALSE; + } + } + + for (int i = 0; i < Npts; i++) { + oprofile[i] = fgaussOD (xprofile[i], par, Npar, NULL); + } + + FREE (par_old); + FREE (wt); +} + +void gaussian_fit (opihi_flt *xprofile, opihi_flt *fprofile, opihi_flt *wprofile, opihi_flt *oprofile, int Npts, opihi_flt *par, int Npar, FitFunc *FUNC) { + + double ochisq = mrqinit (xprofile, fprofile, wprofile, Npts, par, Npar, FUNC, VERBOSE); + double dchisq = ochisq + 2*Npts; + + int Niter; + for (Niter = 0; (Niter < 20) && ((dchisq > 0.1*(Npts - Npar)) || (dchisq <= 0.0)); Niter++) { + double chisq = mrqmin (xprofile, fprofile, wprofile, Npts, par, Npar, FUNC, VERBOSE); + dchisq = ochisq - chisq; + ochisq = chisq; + if (VERBOSE) gprint (GP_ERR, "dchisq: %f, Ndof: %d\n", dchisq, Npts - Npar); + } + if (VERBOSE) gprint (GP_ERR, "%d iterations\n", Niter); + + if (!QUIET) gprint (GP_LOG, "gauss @ %f (%f) : top %f bottom %f\n", par[0], par[1], par[2], par[3]); + + if (oprofile) { + for (int i = 0; i < Npts; i++) { + oprofile[i] = FUNC (xprofile[i], par, Npar, NULL); + } + } + + mrqfree (Npar); +} Index: trunk/Ohana/src/opihi/mana/deimos_fitobj.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_fitobj.c (revision 41158) +++ trunk/Ohana/src/opihi/mana/deimos_fitobj.c (revision 41159) @@ -1,13 +1,4 @@ # include "data.h" # include "deimos.h" - -# if (0) -int deimos_fitobj (int argc, char **argv) { - OHANA_UNUSED_PARAM(argc); - OHANA_UNUSED_PARAM(argv); - return FALSE; -} - -# else /* @@ -18,19 +9,19 @@ // internal functions to fitobj -static float deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts); -int *deimos_sort_result (DeimosResult *result, int Nresult); -int deimos_get_random_sample (int Nresult); -static opihi_flt *deimos_make_test (opihi_flt *guess, float *sigma, int row, int Nrow, int fullInput); - -int USE_GAUSS_DEV = FALSE; +static float deimos_get_chisq (float *buffer, float *vardata, float *model, int Nx, int Ny, int row, int *Npts); +double deimos_LMM_update (float *data, float *vardata, float *model_ref, float *model_obj, float *model_sky, float *model_bck, float *dObj, float *dSky, float *dBck, int Nx, int Ny, int row, float lambda); +static opihi_flt *deimos_make_test (opihi_flt *guess, float sigma, int row, int Nrow, float *dValue); +opihi_flt *deimos_subvector (opihi_flt *fullVec, int row, int Nrow); int deimos_fitobj (int argc, char **argv) { // input parameters: - // * buffer : 2D image of slit (full frame or cutout?) - // * trace : spline fit of slit central x pos vs y-coord - // * profile : slit window profile (vector) - // * PSF : point-spread function vector (flux normalized, x-dir) - // * stilt : slit tilt response : 2D kernel? + // * buffer : 2D image of slit (full frame or cutout?) + // * slit_trace_red : spline fit of slit central x pos vs y-coord + // * slit_trace_blu : spline fit of slit central x pos vs y-coord + // * psf_trace : spline fit of slit central x pos vs y-coord + // * profile : slit window profile (vector) + // * PSF : point-spread function vector (flux normalized, x-dir) + // * stilt : slit tilt response : 2D kernel? // in-out parameters: @@ -41,17 +32,48 @@ int N; - Spline *trace = NULL; - Vector *profile = NULL; - - Vector *obj = NULL; - Vector *sky = NULL; - Vector *bck = NULL; - - Vector *psf = NULL; - - Buffer *buffer = NULL; - - float stilt = 0.0; // angle of the slit + Spline *slit_trace_red = NULL; + Spline *slit_trace_blu = NULL; + Spline *psf_trace = NULL; + + Vector *profile = NULL; + Vector *psf = NULL; + + Vector *obj = NULL; + Vector *sky = NULL; + Vector *bck = NULL; + + Buffer *buffer = NULL; + + float stilt = 0.0; // angle of the slit + ohana_gaussdev_init (); + + int DO_PLUS = FALSE; + if ((N = get_argument (argc, argv, "-plus"))) { + remove_argument (N, &argc, argv); + DO_PLUS = TRUE; + } + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + remove_argument (N, &argc, argv); + VERBOSE = TRUE; + } + + int SKIPFIT = FALSE; + if ((N = get_argument (argc, argv, "-skipfit"))) { + remove_argument (N, &argc, argv); + SKIPFIT = TRUE; + } + + // for a red vs blu spline, we need to specify the split point + // XXX this is REALLY ad-hoc for Deimos. not sure how to make + // this more generic (need to define the ranges somewhere) + int redlimit = 4096; + if ((N = get_argument (argc, argv, "-redlimit"))) { + remove_argument (N, &argc, argv); + redlimit = atoi (argv[N]); + remove_argument (N, &argc, argv); + } // a crude noise model valid for both the guess vectors and the @@ -75,11 +97,10 @@ remove_argument (N, &argc, argv); } - float scale = 0.75; - if ((N = get_argument (argc, argv, "-scale"))) { - remove_argument (N, &argc, argv); - scale = atof(argv[N]); - remove_argument (N, &argc, argv); - } - + float lambda = 1.0; + if ((N = get_argument (argc, argv, "-lambda"))) { + remove_argument (N, &argc, argv); + lambda = atof(argv[N]); + remove_argument (N, &argc, argv); + } int Niter = 25; @@ -89,10 +110,4 @@ remove_argument (N, &argc, argv); } - int Ntry = 5; - if ((N = get_argument (argc, argv, "-try"))) { - remove_argument (N, &argc, argv); - Ntry = atoi(argv[N]); - remove_argument (N, &argc, argv); - } int Nrow = 5; if ((N = get_argument (argc, argv, "-row"))) { @@ -102,26 +117,19 @@ } - USE_GAUSS_DEV = FALSE; - if ((N = get_argument (argc, argv, "-gaussdev"))) { - remove_argument (N, &argc, argv); - USE_GAUSS_DEV = TRUE; - } - - int SAVE_ALL_VECTORS = FALSE; - if ((N = get_argument (argc, argv, "-save-all-vectors"))) { - remove_argument (N, &argc, argv); - SAVE_ALL_VECTORS = TRUE; - } - - int SAVE_MIN_VECTORS = FALSE; - if ((N = get_argument (argc, argv, "-save-min-vectors"))) { - remove_argument (N, &argc, argv); - SAVE_MIN_VECTORS = TRUE; - } - // Input parameters: - if ((N = get_argument (argc, argv, "-trace"))) { - remove_argument (N, &argc, argv); - if ((trace = FindSpline (argv[N])) == NULL) return (FALSE); + if ((N = get_argument (argc, argv, "-slit-trace-red"))) { + remove_argument (N, &argc, argv); + if ((slit_trace_red = FindSpline (argv[N])) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-slit-trace-blu"))) { + remove_argument (N, &argc, argv); + if ((slit_trace_blu = FindSpline (argv[N])) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + + if ((N = get_argument (argc, argv, "-psf-trace"))) { + remove_argument (N, &argc, argv); + if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE); remove_argument (N, &argc, argv); } else { goto usage; } @@ -158,4 +166,20 @@ remove_argument (N, &argc, argv); } else { goto usage; } + + // save the model + Buffer *model = NULL; + if ((N = get_argument (argc, argv, "-model"))) { + remove_argument (N, &argc, argv); + if ((model = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } + + // save the model + Buffer *varim = NULL; + if ((N = get_argument (argc, argv, "-variance"))) { + remove_argument (N, &argc, argv); + if ((varim = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } if (argc != 3) goto usage; @@ -179,4 +203,10 @@ float *bufVal = (float *) buffer->matrix.buffer; + if (varim) { + if (Nx != varim[0].matrix.Naxis[0]) { gprint (GP_ERR, "inconsistent sizes for image and varim\n"); return FALSE; } + if (Ny != varim[0].matrix.Naxis[1]) { gprint (GP_ERR, "inconsistent sizes for image and varim\n"); return FALSE; } + } + float *varVal = varim ? (float *) varim->matrix.buffer : NULL; + // obj, sky, bck must be consistent with data if (Ny != obj->Nelements) { @@ -193,240 +223,141 @@ } - // for the cross-dispersion reference pixle, use the user value or set to Nx/2 if < 0 + if (model) { + if (!CreateBuffer (model, Nx, Ny, -32, 1.0, 0.0)) { gprint (GP_ERR, "error generating output model buffer\n"); return FALSE; } + } + + // for the cross-dispersion reference pixel, use the user value or set to Nx/2 if < 0 deimos_set_cross_ref (atoi(argv[2]), Nx); deimos_make_kernel (stilt, Nx); - // the functions below accept an opihi_flt array for object, sky, background - opihi_flt *objVal = obj->elements.Flt; - opihi_flt *skyVal = sky->elements.Flt; - opihi_flt *bckVal = bck->elements.Flt; - - int NRESULT = Niter + Ntry*Niter + 1; - ALLOCATE_PTR (result, DeimosResult, NRESULT); - - Vector **objOut = NULL; - Vector **skyOut = NULL; - Vector **bckOut = NULL; Vector *objMin = NULL; Vector *skyMin = NULL; Vector *bckMin = NULL; - // Save all result vectors for each of obj, sky, bck: - if (SAVE_ALL_VECTORS) { - ALLOCATE (objOut, Vector *, NRESULT); - ALLOCATE (skyOut, Vector *, NRESULT); - ALLOCATE (bckOut, Vector *, NRESULT); - for (int i = 0; i < NRESULT; i++) { - char name[64]; - snprintf (name, 64, "obj_%04d", i); if ((objOut[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (objOut[i], OPIHI_FLT, Ny); - snprintf (name, 64, "sky_%04d", i); if ((skyOut[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (skyOut[i], OPIHI_FLT, Ny); - snprintf (name, 64, "bck_%04d", i); if ((bckOut[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bckOut[i], OPIHI_FLT, Ny); - for (int iy = 0; iy < Ny; iy++) { - objOut[i]->elements.Flt[iy] = 0.0; - skyOut[i]->elements.Flt[iy] = 0.0; - bckOut[i]->elements.Flt[iy] = 0.0; + if ((objMin = SelectVector ("objMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (objMin, OPIHI_FLT, Ny); + if ((skyMin = SelectVector ("skyMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (skyMin, OPIHI_FLT, Ny); + if ((bckMin = SelectVector ("bckMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bckMin, OPIHI_FLT, Ny); + + // the functions below accept an opihi_flt array for object, sky, background + opihi_flt *objVal = objMin->elements.Flt; + opihi_flt *skyVal = skyMin->elements.Flt; + opihi_flt *bckVal = bckMin->elements.Flt; + + // we make a copy of the input guess vectors and minimize them + for (int iy = 0; iy < Ny; iy++) { + objVal[iy] = obj->elements.Flt[iy]; + skyVal[iy] = sky->elements.Flt[iy]; + bckVal[iy] = bck->elements.Flt[iy]; + } + + float noiseVar = SQ(noise); + ALLOCATE_PTR (objNoise, float, Ny); + ALLOCATE_PTR (skyNoise, float, Ny); + ALLOCATE_PTR (bckNoise, float, Ny); + + // generate noise vectors + for (int iy = 0; iy < Ny; iy++) { + objNoise[iy] = nsigma*sqrt(hypot(fabs(objVal[iy]) / gain, noiseVar)); + skyNoise[iy] = nsigma*sqrt(hypot(fabs(skyVal[iy]) / gain, noiseVar)); + bckNoise[iy] = nsigma*sqrt(hypot(fabs(bckVal[iy]) / gain, noiseVar)); + } + + // obj, sky, bck are initial guesses. Also use these values to define a range for the initial guesses. + + // XXX worry about last segment (< Nrow rows) + for (int iter = 0; iter < Niter; iter++) { + + int dRow = Nrow/2; + + // use scaled-noise to define test ranges + // we are going to try each row using a window Nrow wide to measure chisq: + for (int row = dRow; !SKIPFIT && (row < Ny - dRow - 1); row++) { + + // I am using a numerical equivalent of Levenberg-Marquardt, but only for a single + // p_m[i] = (obj,sky,bck) triplet [m is the index of obj,sky,bck; i is the row + // index] at a time. I start with a current best guess set of (obj,sky,bck). I + // need to calculate the model at this location, and then at offset locations p_m + + // dp_m. Note that I do not need to calculate chisq for each of these to generate + // the elements of the LMM equations: + + int row_0 = row - dRow; // start row of the subset test range + + // Nrow = 5, dRow = 2, Ny = 100, last start is Ny - 2 - 1 = 97 + // last row_0 is 95 + + // generate subset vectors for the range **centered** on row + opihi_flt *objCurr = deimos_subvector (objVal, row_0, Nrow); + opihi_flt *skyCurr = deimos_subvector (skyVal, row_0, Nrow); + opihi_flt *bckCurr = deimos_subvector (bckVal, row_0, Nrow); + + // current vector value: + // for deimos_make_model, row is the starting point of the subimage + float *model_ref = deimos_make_model (objCurr, skyCurr, bckCurr, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Nrow, row_0); + + // in the test functions below, dRow is the central pixel of the subset, e.g., obj[row] => objCurr[dRow] + + float dObj = 0.0, dSky = 0.0, dBck = 0.0; + + // delta obj (save dObj) + opihi_flt *objtest = deimos_make_test (objCurr, objNoise[row], dRow, Nrow, &dObj); + float *model_obj = deimos_make_model (objtest, skyCurr, bckCurr, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Nrow, row_0); + + // delta sky (save dSky) + opihi_flt *skytest = deimos_make_test (skyCurr, skyNoise[row], dRow, Nrow, &dSky); + float *model_sky = deimos_make_model (objCurr, skytest, bckCurr, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Nrow, row_0); + + // delta bck (save dBck) + opihi_flt *bcktest = deimos_make_test (bckCurr, bckNoise[row], dRow, Nrow, &dBck); + float *model_bck = deimos_make_model (objCurr, skyCurr, bcktest, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Nrow, row_0); + + // dObj, dSky, dBck carry the delta used derivatives in and the delta applied to the values out + float chisq = deimos_LMM_update (bufVal, varVal, model_ref, model_obj, model_sky, model_bck, &dObj, &dSky, &dBck, Nx, Nrow, row_0, lambda); + + if (FALSE && (row > 10) && (row < 20)) { + fprintf (stderr, "%d %d : %f %f %f\n", iter, row, dObj, dSky, dBck); } - } - } - if (SAVE_MIN_VECTORS) { - if ((objMin = SelectVector ("objMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (objMin, OPIHI_FLT, Ny); - if ((skyMin = SelectVector ("skyMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (skyMin, OPIHI_FLT, Ny); - if ((bckMin = SelectVector ("bckMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bckMin, OPIHI_FLT, Ny); - for (int iy = 0; iy < Ny; iy++) { - objMin->elements.Flt[iy] = 0.0; - skyMin->elements.Flt[iy] = 0.0; - bckMin->elements.Flt[iy] = 0.0; - } - } - Vector *chisqVect; - if ((chisqVect = SelectVector ("chisqOut", ANYVECTOR, TRUE)) == NULL) { return FALSE; } - ResetVector (chisqVect, OPIHI_FLT, NRESULT); - for (int i = 0; i < NRESULT; i++) { chisqVect->elements.Flt[i] = 0.0; } - - // XXX temp hack: save chisq subset vectors for each row pass - int Npass = Ny / Nrow; - ALLOCATE_PTR (chiPass, Vector *, Npass); - for (int pass = 0; pass < Npass; pass++) { - char name[64]; - snprintf (name, 64, "chiPass_%03d", pass); - if ((chiPass[pass] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return FALSE; } ResetVector (chiPass[pass], OPIHI_FLT, NRESULT); - } - - float noiseVar = SQ(noise); - ALLOCATE_PTR (objNoise, float, Nrow); - ALLOCATE_PTR (skyNoise, float, Nrow); - ALLOCATE_PTR (bckNoise, float, Nrow); - - // work on segments of Nrow at a time - // XXX worry about last segment (< Nrow rows) - int pass = 0; - for (int row = 0; row < Ny - Nrow; row += Nrow, pass++) { - - int Nresult = 0; - - // obj, sky, bck are initial guesses. Also use these values to define a range for the initial guesses. - - // If I just use the value in a given bin to define the range, random noise for the - // faint end could have bad results (range ~ 0.0). For that matter, if I am in a - // region with large dynamic range (e.g., on an emission or absorption line), then we - // likely have large errors in how the guess follows the truth. - - /** old concept: - float objRange = 0, skyRange = 0, bckRange = 0; - for (int iy = 0; iy < Nrow; iy++) { - objRange = MAX(objVal[iy + row], objRange); - skyRange = MAX(skyVal[iy + row], skyRange); - bckRange = MAX(bckVal[iy + row], bckRange); - } - **/ - - // new concept: start with 2-sigma range, then ramp down below - for (int iy = 0; iy < Nrow; iy++) { - objNoise[iy] = nsigma*sqrt(hypot(fabs(objVal[iy + row]) / gain, noiseVar)); - skyNoise[iy] = nsigma*sqrt(hypot(fabs(skyVal[iy + row]) / gain, noiseVar)); - bckNoise[iy] = nsigma*sqrt(hypot(fabs(bckVal[iy + row]) / gain, noiseVar)); - } - - // first test is the input guess - if (1) { - opihi_flt *objtest = deimos_make_test (objVal, NULL, row, Nrow, TRUE); - opihi_flt *skytest = deimos_make_test (skyVal, NULL, row, Nrow, TRUE); - opihi_flt *bcktest = deimos_make_test (bckVal, NULL, row, Nrow, TRUE); - - // generate the model based on the test values - // objtest, skytest, bcktest are segments of the full parameter vectors - float *model = deimos_make_model (objtest, skytest, bcktest, psf, profile, trace, Nx, Nrow, row); - - // calculate the current chisq - // XXX need to add in the per-pixel error or variance - int Npts = 0; - float chisq = deimos_get_chisq (bufVal, model, Nx, Nrow, row, &Npts); - fprintf (stderr, "guess: chisq: %f, Npts: %d\n", chisq, Npts); - free (model); - - if (pass < Npass) { - chiPass[pass]->elements.Flt[Nresult] = chisq / (float) Npts; - } - - // save the results - result[Nresult].chisq = chisq / (float) Npts; - result[Nresult].obj = objtest; - result[Nresult].sky = skytest; - result[Nresult].bck = bcktest; - Nresult ++; - } - - // run Niter iterations - for (int iter = 0; iter < Niter; iter ++) { - - // generate a set of test values for obj, sky, backgnd based on the current guess - opihi_flt *objtest = deimos_make_test (objVal, objNoise, row, Nrow, TRUE); - opihi_flt *skytest = deimos_make_test (skyVal, skyNoise, row, Nrow, TRUE); - opihi_flt *bcktest = deimos_make_test (bckVal, bckNoise, row, Nrow, TRUE); - - for (int itmp = 0; FALSE && (itmp < Nrow); itmp++) { - fprintf (stderr, "val: %f %f %f, tst: %f %f %f\n", - objVal[itmp + row], skyVal[itmp + row], bckVal[itmp + row], - objtest[itmp], skytest[itmp], bcktest[itmp]); - } - - // generate the model based on the test values - // objtest, skytest, bcktest are segments of the full parameter vectors - float *model = deimos_make_model (objtest, skytest, bcktest, psf, profile, trace, Nx, Nrow, row); - - // calculate the current chisq - // XXX need to add in the per-pixel error or variance - int Npts = 0; - float chisq = deimos_get_chisq (bufVal, model, Nx, Nrow, row, &Npts); - // fprintf (stderr, "iter: %d, chisq: %f, Npts: %d\n", iter, chisq, Npts); - free (model); - - if (pass < Npass) { - chiPass[pass]->elements.Flt[Nresult] = chisq / (float) Npts; - } - - // save the results - result[Nresult].chisq = chisq / (float) Npts; - result[Nresult].obj = objtest; - result[Nresult].sky = skytest; - result[Nresult].bck = bcktest; - Nresult ++; - } - - for (int try = 0; try < Ntry; try ++) { - // generate an index of results sorted by chisq - int *IDX = deimos_sort_result (result, Nresult); - - int Nnewset = 0; - - // XXX not sure this scaling factor concept is good - for (int iy = 0; iy < Nrow; iy++) { - objNoise[iy] = scale*objNoise[iy]; - skyNoise[iy] = scale*skyNoise[iy]; - bckNoise[iy] = scale*bckNoise[iy]; - } - - // run Niter iterations, selecting guess values from the result set - for (int iter = 0; iter < Niter; iter ++) { - int bin = deimos_get_random_sample (Nresult); - int entry = IDX[bin]; - // fprintf (stderr, "try: %d : %d, %f\n", bin, entry, result[entry].chisq); - - // generate a set of test values for obj, sky, backgnd based on the current guess - opihi_flt *objtest = deimos_make_test (result[entry].obj, objNoise, row, Nrow, FALSE); - opihi_flt *skytest = deimos_make_test (result[entry].sky, skyNoise, row, Nrow, FALSE); - opihi_flt *bcktest = deimos_make_test (result[entry].bck, bckNoise, row, Nrow, FALSE); - - for (int itmp = 0; FALSE && (itmp < Nrow); itmp++) { - fprintf (stderr, "val: %f %f %f, tst: %f %f %f\n", - result[entry].obj[itmp], result[entry].sky[itmp], result[entry].bck[itmp], - objtest[itmp], skytest[itmp], bcktest[itmp]); - } - - // generate the model based on the test values - float *model = deimos_make_model (objtest, skytest, bcktest, psf, profile, trace, Nx, Nrow, row); - - // calculate the current chisq - int Npts = 0; - float chisq = deimos_get_chisq (bufVal, model, Nx, Nrow, row, &Npts); - // fprintf (stderr, "try: %d, iter: %d, chisq: %f, Npts: %d\n", try, iter, chisq, Npts); - free (model); - - if (pass < Npass) { - chiPass[pass]->elements.Flt[Nresult + Nnewset] = chisq / (float) Npts; - } - - // save the results - result[Nresult + Nnewset].chisq = chisq / (float) Npts; - result[Nresult + Nnewset].obj = objtest; - result[Nresult + Nnewset].sky = skytest; - result[Nresult + Nnewset].bck = bcktest; - Nnewset ++; - } - Nresult += Nnewset; - } - - int *IDX = deimos_sort_result (result, Nresult); - - // ENDING: save the best result and info about the allowed chisq region - // int *index = sort_result (result, Nresult); - // copy result[index[0]].obj,sky,bck to obj,sky,bck - for (int i = 0; i < Nresult; i++) { - for (int iy = 0; iy < Nrow; iy++) { - if (SAVE_ALL_VECTORS) { - objOut[i]->elements.Flt[iy + row] = result[IDX[i]].obj[iy]; - skyOut[i]->elements.Flt[iy + row] = result[IDX[i]].sky[iy]; - bckOut[i]->elements.Flt[iy + row] = result[IDX[i]].bck[iy]; + + if (1) { + if (DO_PLUS) { + objVal[row] += dObj; + skyVal[row] += dSky; + bckVal[row] += dBck; + } else { + objVal[row] -= dObj; + skyVal[row] -= dSky; + bckVal[row] -= dBck; } } - chisqVect->elements.Flt[i] += result[IDX[i]].chisq; + + free (objCurr); + free (skyCurr); + free (bckCurr); + + free (model_ref); + + free (model_obj); + free (model_sky); + free (model_bck); + + free (objtest); + free (skytest); + free (bcktest); + + if (VERBOSE) fprintf (stderr, "chisq: %f (%f, %f, %f)\n", chisq, dObj, dSky, dBck); } - if (SAVE_MIN_VECTORS) { - for (int iy = 0; iy < Nrow; iy++) { - objMin->elements.Flt[iy + row] = result[IDX[0]].obj[iy]; - skyMin->elements.Flt[iy + row] = result[IDX[0]].sky[iy]; - bckMin->elements.Flt[iy + row] = result[IDX[0]].bck[iy]; - } + + float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Ny, 0); + + int Npts = 0; + float chisq_full = deimos_get_chisq (bufVal, varVal, model_full, Nx, Ny, 0, &Npts); + + if (model) { + free (model[0].matrix.buffer); + model[0].matrix.buffer = (char *) model_full; + } else { + free (model_full); } + + fprintf (stderr, "** Full chisq: %f (%d) : %f)\n", chisq_full, Npts, chisq_full / (1.0 * Npts)); } @@ -434,5 +365,6 @@ usage: - gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n"); + gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -slit-trace-red (spline) -slit-trace-blu (spline) -psf-trace (spline) -profile vector -psf vector -stilt (angle)\n"); + gprint (GP_ERR, "OPTIONS: [-model image] [-variance image] [-gain value] [-noise value] [-nsigma value] [-lambda value] [-iter N] [-row N] [-redlimit Npixel] [-skipfit] [-v] [-plus]\n"); return FALSE; } @@ -440,38 +372,120 @@ /****************** fitobj Support Functions *******************/ -static opihi_flt *deimos_make_test (opihi_flt *guess, float *sigma, int row, int Nrow, int fullInput) { - - // we have a vector, guess, with some length >= row + Nrow. - // generate a set of random values in the vicinty of guess[row] to guess[row+Nrow] - // return vector is of length Nrow +// we generate a new vector with a single element (row) modified based on sigma +opihi_flt *deimos_subvector (opihi_flt *fullVec, int row, int Nrow) { + + ALLOCATE_PTR (value, opihi_flt, Nrow); + + for (int i = 0; i < Nrow; i++) { + value[i] = fullVec[i + row]; + } + + return value; +} + + +// we generate a new vector with a single element (row) modified based on sigma +static opihi_flt *deimos_make_test (opihi_flt *guess, float sigma, int row, int Nrow, float *dValue) { // if fullInput is TRUE, the guess vector runs from 0 to Ny (full wavelength range) // if fullInput is FALSE, the guess vector runs from row to row + Nrow (subset range) - // sigma is always only Nrow long + // sigma is valid for the element at row ALLOCATE_PTR (value, opihi_flt, Nrow); - int start_pix = fullInput ? row : 0; - int stop_pix = fullInput ? row + Nrow : Nrow; - - int iv = 0; - for (int i = start_pix; i < stop_pix; i++, iv++) { - if (!sigma) { - value[iv] = guess[i]; - continue; - } - if (USE_GAUSS_DEV) { - value[iv] = ohana_gaussdev_rnd(guess[i], sigma[iv]); - continue; - } - - // uniform distribution: - value[iv] = sigma[iv]*(drand48() - 0.5) + guess[i]; - } + for (int i = 0; i < Nrow; i++) { + value[i] = guess[i]; + } + + value[row] = guess[row] + 0.1*sigma; + + *dValue = 0.1*sigma; return value; } -static float deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts) { +double deimos_LMM_update (float *data, float *vardata, float *model_ref, + float *model_obj, float *model_sky, float *model_bck, + float *dObj, float *dSky, float *dBck, int Nx, int Ny, int row, float lambda) { + + // we are going to loop over the images calculating the following value: + + // W = 1 / sigma^2, leave as 1.0 for now + // chisq = sum W * (data - model_ref)^2 + // dF_X = sum W * (data - model_ref) * (model_ref - model_X) / dX + // d2F_XY = sum W * (model_ref - model_X) * (model_ref - model_Y) / (dX * dY) + + double chisq = 0.0; + double dF_obj = 0.0, dF_sky = 0.0, dF_bck = 0.0; + double d2F_obj_obj = 0.0, d2F_obj_sky = 0.0, d2F_obj_bck = 0.0, d2F_sky_sky = 0.0, d2F_sky_bck = 0.0, d2F_bck_bck = 0.0; + + for (int iy = 0; iy < Ny; iy++) { + for (int ix = 0; ix < Nx; ix++) { + int pix = ix + iy*Nx; + int pix_data = ix + (iy + row)*Nx; + + // if any of the pixel values are NAN, skip the point: + if (!isfinite(data[pix_data])) continue; + if (!isfinite(model_ref[pix])) continue; + if (!isfinite(model_obj[pix])) continue; + if (!isfinite(model_sky[pix])) continue; + if (!isfinite(model_bck[pix])) continue; + + if (vardata && !isfinite(vardata[pix_data])) continue; + float Var = vardata ? vardata[pix_data] : 1.0; + float W = (fabs(Var) < 0.01) ? 100.0 : 1.0 / Var; + + float df = W * (data[pix_data] - model_ref[pix]); + + chisq += SQ(df); + + float df_obj = (model_ref[pix] - model_obj[pix]) / *dObj; + float df_sky = (model_ref[pix] - model_sky[pix]) / *dSky; + float df_bck = (model_ref[pix] - model_bck[pix]) / *dBck; + + dF_obj += df*df_obj; + dF_sky += df*df_sky; + dF_bck += df*df_bck; + + d2F_obj_obj += W * SQ(df_obj); + d2F_sky_sky += W * SQ(df_sky); + d2F_bck_bck += W * SQ(df_bck); + + d2F_obj_sky += W * df_obj * df_sky; + d2F_obj_bck += W * df_obj * df_bck; + d2F_sky_bck += W * df_sky * df_bck; + } + } + + // these are the elements of Ax = B + ALLOCATE_PTR (A, double *, 3); + ALLOCATE_PTR (B, double *, 3); + for (int i = 0; i < 3; i++) { + ALLOCATE (A[i], double, 3); + ALLOCATE (B[i], double, 1); + } + + A[0][0] = d2F_obj_obj * (1.0 + lambda); + A[1][1] = d2F_sky_sky * (1.0 + lambda); + A[2][2] = d2F_bck_bck * (1.0 + lambda); + + A[0][1] = A[1][0] = d2F_obj_sky; + A[1][2] = A[2][1] = d2F_sky_bck; + A[0][2] = A[2][0] = d2F_obj_bck; + + B[0][0] = dF_obj; + B[1][0] = dF_sky; + B[2][0] = dF_bck; + + dgaussjordan (A, B, 3, 1); + + *dObj = B[0][0]; + *dSky = B[1][0]; + *dBck = B[2][0]; + + return chisq; +} + +static float deimos_get_chisq (float *buffer, float *vardata, float *model, int Nx, int Ny, int row, int *Npts) { int npts = 0; @@ -487,5 +501,12 @@ if (!isfinite(buffer[pix_buffer])) continue; if (!isfinite(model[pix_model])) continue; - chisq += SQ(buffer[pix_buffer] - model[pix_model]); + + if (vardata && !isfinite(vardata[pix_buffer])) continue; + float Var = vardata ? vardata[pix_buffer] : 1.0; + float W = (fabs(Var) < 0.01) ? 100.0 : 1.0 / Var; + + float dchisq = SQ(buffer[pix_buffer] - model[pix_model]) * W; + if (!isfinite(dchisq)) continue; + chisq += dchisq; npts ++; } @@ -494,42 +515,2 @@ return chisq; } - -int *deimos_sort_result (DeimosResult *result, int Nresult) { - - ALLOCATE_PTR (IDX, int, Nresult); - ALLOCATE_PTR (chisq, float, Nresult); - - for (int i = 0; i < Nresult; i++) { - IDX[i] = i; - chisq[i] = result[i].chisq; - } - - sort_float_index (chisq, IDX, Nresult); - free (chisq); - - return IDX; -} - -static float DEIMOS_A = 2.0; // value for random sample - -int deimos_get_random_sample (int Nresult) { - - // choose a bin from the range 0 - Nresult, with front-loaded weighting. - - // bin = Nresult * (exp(q) - exp(-A)) / (exp(A) - exp(-A)) - // q = 2A(x - 1/2) : random number between -A and +A - // x : random number between 0 and 1 - - // f(q;A) = (exp(q) - exp(-A)) / (exp(A) - exp(-A)) - // f(q;A) ranges from 0 to +1 as x ranges from 0 to +1 - - float x = drand48(); - float q = 2*DEIMOS_A*(x - 0.5); - float f = (exp(q) - exp(-DEIMOS_A)) / (exp(DEIMOS_A) - exp(-DEIMOS_A)); - // XXX optimization : precalculate and save the elements of this function - - int bin = floor(Nresult * f); - - return bin; -} -# endif Index: trunk/Ohana/src/opihi/mana/deimos_fitobj.v0.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_fitobj.v0.c (revision 41159) +++ trunk/Ohana/src/opihi/mana/deimos_fitobj.v0.c (revision 41159) @@ -0,0 +1,538 @@ +# include "data.h" +# include "deimos.h" + +# if (0) +int deimos_fitobj (int argc, char **argv) { + OHANA_UNUSED_PARAM(argc); + OHANA_UNUSED_PARAM(argv); + return FALSE; +} + +# else + +/* + this is starting to work OK. some improvements to make + * use Gaussdev to sample + * do not scale down range (or user-set scale-down) + */ + +// internal functions to fitobj +static float deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts); +int *deimos_sort_result (DeimosResult *result, int Nresult); +int deimos_get_random_sample (int Nresult); +static opihi_flt *deimos_make_test (opihi_flt *guess, float *sigma, int row, int Nrow, int fullInput); + +int USE_GAUSS_DEV = FALSE; + +int deimos_fitobj (int argc, char **argv) { + + // input parameters: + // * buffer : 2D image of slit (full frame or cutout?) + // * trace : spline fit of slit central x pos vs y-coord + // * profile : slit window profile (vector) + // * PSF : point-spread function vector (flux normalized, x-dir) + // * stilt : slit tilt response : 2D kernel? + + // in-out parameters: + // * obj : vector of object flux vs y-coord (starting guess and result) + // * sky : vector of local sky signal vs y-coord + // * bck : vector of extra-slit background flux vs y-coord + + int N; + + Spline *trace = NULL; + Vector *profile = NULL; + + Vector *obj = NULL; + Vector *sky = NULL; + Vector *bck = NULL; + + Vector *psf = NULL; + + Buffer *buffer = NULL; + + float stilt = 0.0; // angle of the slit + ohana_gaussdev_init (); + + fprintf (stderr, "use deimos fitalt instead\n"); + return FALSE; + + // a crude noise model valid for both the guess vectors and the + // data buffer: + float gain = 1.0; + if ((N = get_argument (argc, argv, "-gain"))) { + remove_argument (N, &argc, argv); + gain = atof(argv[N]); + remove_argument (N, &argc, argv); + } + float noise = 0.0; + if ((N = get_argument (argc, argv, "-noise"))) { + remove_argument (N, &argc, argv); + noise = atof(argv[N]); + remove_argument (N, &argc, argv); + } + float nsigma = 2.0; + if ((N = get_argument (argc, argv, "-nsigma"))) { + remove_argument (N, &argc, argv); + nsigma = atof(argv[N]); + remove_argument (N, &argc, argv); + } + float scale = 0.75; + if ((N = get_argument (argc, argv, "-scale"))) { + remove_argument (N, &argc, argv); + scale = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + + int Niter = 25; + if ((N = get_argument (argc, argv, "-iter"))) { + remove_argument (N, &argc, argv); + Niter = atoi(argv[N]); + remove_argument (N, &argc, argv); + } + int Ntry = 5; + if ((N = get_argument (argc, argv, "-try"))) { + remove_argument (N, &argc, argv); + Ntry = atoi(argv[N]); + remove_argument (N, &argc, argv); + } + int Nrow = 5; + if ((N = get_argument (argc, argv, "-row"))) { + remove_argument (N, &argc, argv); + Nrow = atoi(argv[N]); + remove_argument (N, &argc, argv); + } + + USE_GAUSS_DEV = FALSE; + if ((N = get_argument (argc, argv, "-gaussdev"))) { + remove_argument (N, &argc, argv); + USE_GAUSS_DEV = TRUE; + } + + int SAVE_ALL_VECTORS = FALSE; + if ((N = get_argument (argc, argv, "-save-all-vectors"))) { + remove_argument (N, &argc, argv); + SAVE_ALL_VECTORS = TRUE; + } + + int SAVE_MIN_VECTORS = FALSE; + if ((N = get_argument (argc, argv, "-save-min-vectors"))) { + remove_argument (N, &argc, argv); + SAVE_MIN_VECTORS = TRUE; + } + + // Input parameters: + if ((N = get_argument (argc, argv, "-trace"))) { + remove_argument (N, &argc, argv); + if ((trace = FindSpline (argv[N])) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-profile"))) { + remove_argument (N, &argc, argv); + if ((profile = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-psf"))) { + remove_argument (N, &argc, argv); + if ((psf = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-stilt"))) { + remove_argument (N, &argc, argv); + stilt = atof (argv[N]); + remove_argument (N, &argc, argv); + } else { goto usage; } + + // In-Out parameters: + if ((N = get_argument (argc, argv, "-object"))) { + remove_argument (N, &argc, argv); + if ((obj = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-sky"))) { + remove_argument (N, &argc, argv); + if ((sky = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-backgnd"))) { + remove_argument (N, &argc, argv); + if ((bck = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + + if (argc != 3) goto usage; + + // buffer is the full image, Xref is reference coordinate of the profile in the buffer + if ((buffer = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + + // profile and PSF are defined with central reference pixel mapped to Xref (+ trace) + if (profile->Nelements % 2 == 0) { + gprint (GP_ERR, "slit profile vector must have an odd number of pixels\n"); + return FALSE; + } + if (psf->Nelements % 2 == 0) { + gprint (GP_ERR, "PSF vector must have an odd number of pixels\n"); + return FALSE; + } + + // observed buffer size + int Nx = buffer[0].matrix.Naxis[0]; + int Ny = buffer[0].matrix.Naxis[1]; + float *bufVal = (float *) buffer->matrix.buffer; + + // obj, sky, bck must be consistent with data + if (Ny != obj->Nelements) { + gprint (GP_ERR, "inconsistent wavelength scales (object)\n"); + return FALSE; + } + if (Ny != sky->Nelements) { + gprint (GP_ERR, "inconsistent wavelength scales (sky)\n"); + return FALSE; + } + if (Ny != bck->Nelements) { + gprint (GP_ERR, "inconsistent wavelength scales (backgnd)\n"); + return FALSE; + } + + // for the cross-dispersion reference pixle, use the user value or set to Nx/2 if < 0 + deimos_set_cross_ref (atoi(argv[2]), Nx); + deimos_make_kernel (stilt, Nx); + + // the functions below accept an opihi_flt array for object, sky, background + opihi_flt *objVal = obj->elements.Flt; + opihi_flt *skyVal = sky->elements.Flt; + opihi_flt *bckVal = bck->elements.Flt; + + int NRESULT = Niter + Ntry*Niter + 1; + ALLOCATE_PTR (result, DeimosResult, NRESULT); + + Vector **objOut = NULL; + Vector **skyOut = NULL; + Vector **bckOut = NULL; + Vector *objMin = NULL; + Vector *skyMin = NULL; + Vector *bckMin = NULL; + + // Save all result vectors for each of obj, sky, bck: + if (SAVE_ALL_VECTORS) { + ALLOCATE (objOut, Vector *, NRESULT); + ALLOCATE (skyOut, Vector *, NRESULT); + ALLOCATE (bckOut, Vector *, NRESULT); + for (int i = 0; i < NRESULT; i++) { + char name[64]; + snprintf (name, 64, "obj_%04d", i); if ((objOut[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (objOut[i], OPIHI_FLT, Ny); + snprintf (name, 64, "sky_%04d", i); if ((skyOut[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (skyOut[i], OPIHI_FLT, Ny); + snprintf (name, 64, "bck_%04d", i); if ((bckOut[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bckOut[i], OPIHI_FLT, Ny); + for (int iy = 0; iy < Ny; iy++) { + objOut[i]->elements.Flt[iy] = 0.0; + skyOut[i]->elements.Flt[iy] = 0.0; + bckOut[i]->elements.Flt[iy] = 0.0; + } + } + } + if (SAVE_MIN_VECTORS) { + if ((objMin = SelectVector ("objMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (objMin, OPIHI_FLT, Ny); + if ((skyMin = SelectVector ("skyMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (skyMin, OPIHI_FLT, Ny); + if ((bckMin = SelectVector ("bckMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bckMin, OPIHI_FLT, Ny); + for (int iy = 0; iy < Ny; iy++) { + objMin->elements.Flt[iy] = 0.0; + skyMin->elements.Flt[iy] = 0.0; + bckMin->elements.Flt[iy] = 0.0; + } + } + Vector *chisqVect; + if ((chisqVect = SelectVector ("chisqOut", ANYVECTOR, TRUE)) == NULL) { return FALSE; } + ResetVector (chisqVect, OPIHI_FLT, NRESULT); + for (int i = 0; i < NRESULT; i++) { chisqVect->elements.Flt[i] = 0.0; } + + // XXX temp hack: save chisq subset vectors for each row pass + int Npass = Ny / Nrow; + ALLOCATE_PTR (chiPass, Vector *, Npass); + for (int pass = 0; pass < Npass; pass++) { + char name[64]; + snprintf (name, 64, "chiPass_%03d", pass); + if ((chiPass[pass] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return FALSE; } ResetVector (chiPass[pass], OPIHI_FLT, NRESULT); + } + + float noiseVar = SQ(noise); + ALLOCATE_PTR (objNoise, float, Nrow); + ALLOCATE_PTR (skyNoise, float, Nrow); + ALLOCATE_PTR (bckNoise, float, Nrow); + + // work on segments of Nrow at a time + // XXX worry about last segment (< Nrow rows) + int pass = 0; + for (int row = 0; row < Ny - Nrow; row += Nrow, pass++) { + + int Nresult = 0; + + // obj, sky, bck are initial guesses. Also use these values to define a range for the initial guesses. + + // If I just use the value in a given bin to define the range, random noise for the + // faint end could have bad results (range ~ 0.0). For that matter, if I am in a + // region with large dynamic range (e.g., on an emission or absorption line), then we + // likely have large errors in how the guess follows the truth. + + /** old concept: + float objRange = 0, skyRange = 0, bckRange = 0; + for (int iy = 0; iy < Nrow; iy++) { + objRange = MAX(objVal[iy + row], objRange); + skyRange = MAX(skyVal[iy + row], skyRange); + bckRange = MAX(bckVal[iy + row], bckRange); + } + **/ + + // new concept: start with 2-sigma range, then ramp down below + for (int iy = 0; iy < Nrow; iy++) { + objNoise[iy] = nsigma*sqrt(hypot(fabs(objVal[iy + row]) / gain, noiseVar)); + skyNoise[iy] = nsigma*sqrt(hypot(fabs(skyVal[iy + row]) / gain, noiseVar)); + bckNoise[iy] = nsigma*sqrt(hypot(fabs(bckVal[iy + row]) / gain, noiseVar)); + } + + // first test is the input guess + if (1) { + opihi_flt *objtest = deimos_make_test (objVal, NULL, row, Nrow, TRUE); + opihi_flt *skytest = deimos_make_test (skyVal, NULL, row, Nrow, TRUE); + opihi_flt *bcktest = deimos_make_test (bckVal, NULL, row, Nrow, TRUE); + + // generate the model based on the test values + // objtest, skytest, bcktest are segments of the full parameter vectors + float *model = deimos_make_model (objtest, skytest, bcktest, psf, profile, trace, NULL, Nx, Nrow, row); + + // calculate the current chisq + // XXX need to add in the per-pixel error or variance + int Npts = 0; + float chisq = deimos_get_chisq (bufVal, model, Nx, Nrow, row, &Npts); + fprintf (stderr, "guess: chisq: %f, Npts: %d\n", chisq, Npts); + free (model); + + if (pass < Npass) { + chiPass[pass]->elements.Flt[Nresult] = chisq / (float) Npts; + } + + // save the results + result[Nresult].chisq = chisq / (float) Npts; + result[Nresult].obj = objtest; + result[Nresult].sky = skytest; + result[Nresult].bck = bcktest; + Nresult ++; + } + + // run Niter iterations + for (int iter = 0; iter < Niter; iter ++) { + + // generate a set of test values for obj, sky, backgnd based on the current guess + opihi_flt *objtest = deimos_make_test (objVal, objNoise, row, Nrow, TRUE); + opihi_flt *skytest = deimos_make_test (skyVal, skyNoise, row, Nrow, TRUE); + opihi_flt *bcktest = deimos_make_test (bckVal, bckNoise, row, Nrow, TRUE); + + for (int itmp = 0; FALSE && (itmp < Nrow); itmp++) { + fprintf (stderr, "val: %f %f %f, tst: %f %f %f\n", + objVal[itmp + row], skyVal[itmp + row], bckVal[itmp + row], + objtest[itmp], skytest[itmp], bcktest[itmp]); + } + + // generate the model based on the test values + // objtest, skytest, bcktest are segments of the full parameter vectors + float *model = deimos_make_model (objtest, skytest, bcktest, psf, profile, trace, NULL, Nx, Nrow, row); + + // calculate the current chisq + // XXX need to add in the per-pixel error or variance + int Npts = 0; + float chisq = deimos_get_chisq (bufVal, model, Nx, Nrow, row, &Npts); + // fprintf (stderr, "iter: %d, chisq: %f, Npts: %d\n", iter, chisq, Npts); + free (model); + + if (pass < Npass) { + chiPass[pass]->elements.Flt[Nresult] = chisq / (float) Npts; + } + + // save the results + result[Nresult].chisq = chisq / (float) Npts; + result[Nresult].obj = objtest; + result[Nresult].sky = skytest; + result[Nresult].bck = bcktest; + Nresult ++; + } + + for (int try = 0; try < Ntry; try ++) { + // generate an index of results sorted by chisq + int *IDX = deimos_sort_result (result, Nresult); + + int Nnewset = 0; + + // XXX not sure this scaling factor concept is good + for (int iy = 0; iy < Nrow; iy++) { + objNoise[iy] = scale*objNoise[iy]; + skyNoise[iy] = scale*skyNoise[iy]; + bckNoise[iy] = scale*bckNoise[iy]; + } + + // run Niter iterations, selecting guess values from the result set + for (int iter = 0; iter < Niter; iter ++) { + int bin = deimos_get_random_sample (Nresult); + int entry = IDX[bin]; + // fprintf (stderr, "try: %d : %d, %f\n", bin, entry, result[entry].chisq); + + // generate a set of test values for obj, sky, backgnd based on the current guess + opihi_flt *objtest = deimos_make_test (result[entry].obj, objNoise, row, Nrow, FALSE); + opihi_flt *skytest = deimos_make_test (result[entry].sky, skyNoise, row, Nrow, FALSE); + opihi_flt *bcktest = deimos_make_test (result[entry].bck, bckNoise, row, Nrow, FALSE); + + for (int itmp = 0; FALSE && (itmp < Nrow); itmp++) { + fprintf (stderr, "val: %f %f %f, tst: %f %f %f\n", + result[entry].obj[itmp], result[entry].sky[itmp], result[entry].bck[itmp], + objtest[itmp], skytest[itmp], bcktest[itmp]); + } + + // generate the model based on the test values + float *model = deimos_make_model (objtest, skytest, bcktest, psf, profile, trace, NULL, Nx, Nrow, row); + + // calculate the current chisq + int Npts = 0; + float chisq = deimos_get_chisq (bufVal, model, Nx, Nrow, row, &Npts); + // fprintf (stderr, "try: %d, iter: %d, chisq: %f, Npts: %d\n", try, iter, chisq, Npts); + free (model); + + if (pass < Npass) { + chiPass[pass]->elements.Flt[Nresult + Nnewset] = chisq / (float) Npts; + } + + // save the results + result[Nresult + Nnewset].chisq = chisq / (float) Npts; + result[Nresult + Nnewset].obj = objtest; + result[Nresult + Nnewset].sky = skytest; + result[Nresult + Nnewset].bck = bcktest; + Nnewset ++; + } + Nresult += Nnewset; + } + + int *IDX = deimos_sort_result (result, Nresult); + + // ENDING: save the best result and info about the allowed chisq region + // int *index = sort_result (result, Nresult); + // copy result[index[0]].obj,sky,bck to obj,sky,bck + for (int i = 0; i < Nresult; i++) { + for (int iy = 0; iy < Nrow; iy++) { + if (SAVE_ALL_VECTORS) { + objOut[i]->elements.Flt[iy + row] = result[IDX[i]].obj[iy]; + skyOut[i]->elements.Flt[iy + row] = result[IDX[i]].sky[iy]; + bckOut[i]->elements.Flt[iy + row] = result[IDX[i]].bck[iy]; + } + } + chisqVect->elements.Flt[i] += result[IDX[i]].chisq; + } + if (SAVE_MIN_VECTORS) { + for (int iy = 0; iy < Nrow; iy++) { + objMin->elements.Flt[iy + row] = result[IDX[0]].obj[iy]; + skyMin->elements.Flt[iy + row] = result[IDX[0]].sky[iy]; + bckMin->elements.Flt[iy + row] = result[IDX[0]].bck[iy]; + } + } + } + + return TRUE; + + usage: + gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n"); + return FALSE; +} + +/****************** fitobj Support Functions *******************/ + +static opihi_flt *deimos_make_test (opihi_flt *guess, float *sigma, int row, int Nrow, int fullInput) { + + // we have a vector, guess, with some length >= row + Nrow. + // generate a set of random values in the vicinty of guess[row] to guess[row+Nrow] + // return vector is of length Nrow + + // if fullInput is TRUE, the guess vector runs from 0 to Ny (full wavelength range) + // if fullInput is FALSE, the guess vector runs from row to row + Nrow (subset range) + // sigma is always only Nrow long + + ALLOCATE_PTR (value, opihi_flt, Nrow); + + int start_pix = fullInput ? row : 0; + int stop_pix = fullInput ? row + Nrow : Nrow; + + int iv = 0; + for (int i = start_pix; i < stop_pix; i++, iv++) { + if (!sigma) { + value[iv] = guess[i]; + continue; + } + if (USE_GAUSS_DEV) { + value[iv] = ohana_gaussdev_rnd(guess[i], sigma[iv]); + continue; + } + + // uniform distribution: + value[iv] = sigma[iv]*(drand48() - 0.5) + guess[i]; + } + return value; +} + + +static float deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts) { + + int npts = 0; + float chisq = 0; + + for (int iy = 0; iy < Ny; iy++) { + + for (int ix = 0; ix < Nx; ix++) { + + int pix_buffer = ix + (iy + row)*Nx; + int pix_model = ix + iy*Nx; + + if (!isfinite(buffer[pix_buffer])) continue; + if (!isfinite(model[pix_model])) continue; + chisq += SQ(buffer[pix_buffer] - model[pix_model]); + npts ++; + } + } + *Npts = npts; + return chisq; +} + +int *deimos_sort_result (DeimosResult *result, int Nresult) { + + ALLOCATE_PTR (IDX, int, Nresult); + ALLOCATE_PTR (chisq, float, Nresult); + + for (int i = 0; i < Nresult; i++) { + IDX[i] = i; + chisq[i] = result[i].chisq; + } + + sort_float_index (chisq, IDX, Nresult); + free (chisq); + + return IDX; +} + +static float DEIMOS_A = 2.0; // value for random sample + +int deimos_get_random_sample (int Nresult) { + + // choose a bin from the range 0 - Nresult, with front-loaded weighting. + + // bin = Nresult * (exp(q) - exp(-A)) / (exp(A) - exp(-A)) + // q = 2A(x - 1/2) : random number between -A and +A + // x : random number between 0 and 1 + + // f(q;A) = (exp(q) - exp(-A)) / (exp(A) - exp(-A)) + // f(q;A) ranges from 0 to +1 as x ranges from 0 to +1 + + float x = drand48(); + float q = 2*DEIMOS_A*(x - 0.5); + float f = (exp(q) - exp(-DEIMOS_A)) / (exp(DEIMOS_A) - exp(-DEIMOS_A)); + // XXX optimization : precalculate and save the elements of this function + + int bin = floor(Nresult * f); + + return bin; +} +# endif Index: trunk/Ohana/src/opihi/mana/deimos_fitprofile.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_fitprofile.c (revision 41159) +++ trunk/Ohana/src/opihi/mana/deimos_fitprofile.c (revision 41159) @@ -0,0 +1,299 @@ +# include "data.h" + +typedef opihi_flt FitFunc (opihi_flt, opihi_flt *, int, opihi_flt *); +typedef enum {MODE_UPPER, MODE_LOWER} SigmoidMode; + +void sigmoid_fit (opihi_flt *xprofile, opihi_flt *fprofile, opihi_flt *wprofile, opihi_flt *oprofile, int Npts, opihi_flt *par, int Npar, FitFunc *FUNC); +void sigmoid_fit_irls (opihi_flt *xprofile, opihi_flt *fprofile, opihi_flt *wprofile, opihi_flt *oprofile, int NptsAll, opihi_flt *par, int Npar, SigmoidMode mode, double limit); + +opihi_flt fsigmoidDeimos (opihi_flt, opihi_flt *, int, opihi_flt *); +opihi_flt rsigmoidDeimos (opihi_flt, opihi_flt *, int, opihi_flt *); + +double deimos_weight_cauchy (double x); + +// use these globals to carry into sigmoid_fit +static int VERBOSE = FALSE; +static int QUIET = FALSE; +static int MaxIterations = 10; + +int deimos_fitprofile (int argc, char **argv) { + + // fit simple slit window profile model to vector + // I am fitting a sigmoid at both edges + + // input parameters: + // * xprofile : extracted slit window coordinates (vector) + // * xprofile : extracted slit window profile (vector) + // * Xs, Xe : scale coordinates of window boundaries + + // output parameters: + // * sigmoid_left_center, etc : coordinates of window boundaries + + int N; + Vector *xprofile = NULL; + Vector *fprofile = NULL; + Vector *wprofile = NULL; + + QUIET = FALSE; + if ((N = get_argument (argc, argv, "-q"))) { + QUIET = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-quiet"))) { + QUIET = TRUE; + remove_argument (N, &argc, argv); + } + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + MaxIterations = 10; + if ((N = get_argument (argc, argv, "-max-iterations"))) { + remove_argument (N, &argc, argv); + MaxIterations = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 6) goto usage; + if ((xprofile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((fprofile = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((wprofile = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if (xprofile->Nelements != fprofile->Nelements) goto usage; + if (xprofile->Nelements != wprofile->Nelements) goto usage; + + // supplied guess for boundaries + double Xs = atof(argv[4]); + double Xe = atof(argv[5]); + + // output vectors of fits + Vector *Lprofile = SelectVector ("sigmoid_left", ANYVECTOR, TRUE); + if (Lprofile == NULL) return (FALSE); + MatchVector (Lprofile, xprofile, OPIHI_FLT); + Vector *Rprofile = SelectVector ("sigmoid_right", ANYVECTOR, TRUE); + if (Rprofile == NULL) return (FALSE); + MatchVector (Rprofile, xprofile, OPIHI_FLT); + + // first, use supplied Xs & Xe guess to make a guess for BckL, Sky, BckR + + float Sky_V = 0.0, BckL_V = 0.0, BckR_V = 0.0; + float Sky_N = 0.0, BckL_N = 0.0, BckR_N = 0.0; + for (int i = 0; i < xprofile->Nelements; i++) { + if (!isfinite(wprofile->elements.Flt[i])) continue; + if (!isfinite(fprofile->elements.Flt[i])) continue; + float wt = 1 / SQ(wprofile->elements.Flt[i]); + if (xprofile->elements.Flt[i] < Xs) { + BckL_V += fprofile->elements.Flt[i]*wt; + BckL_N += wt; + continue; + } + if (xprofile->elements.Flt[i] > Xe) { + BckR_V += fprofile->elements.Flt[i]*wt; + BckR_N += wt; + continue; + } + Sky_V += fprofile->elements.Flt[i]*wt; + Sky_N += wt; + } + float Sky = Sky_V / Sky_N; + float BckL = BckL_V / BckL_N; + float BckR = BckR_V / BckR_N; + + // temporary storage vectors + opihi_flt par[4]; + + float Xmid = 0.5*(Xs + Xe); + + int Npar = 4; + par[0] = Xs; + par[1] = 1; + par[2] = Sky; + par[3] = BckL; + sigmoid_fit_irls (xprofile->elements.Flt, fprofile->elements.Flt, wprofile->elements.Flt, Lprofile->elements.Flt, xprofile->Nelements, par, Npar, MODE_LOWER, Xmid); + set_variable ("sigmoid_left_center", par[0]); + set_variable ("sigmoid_left_sigma", par[1]); + set_variable ("sigmoid_left_top", par[2]); + set_variable ("sigmoid_left_bottom", par[3]); + + par[0] = Xe; + par[1] = 1; + par[2] = Sky; + par[3] = BckR; + sigmoid_fit_irls (xprofile->elements.Flt, fprofile->elements.Flt, wprofile->elements.Flt, Rprofile->elements.Flt, xprofile->Nelements, par, Npar, MODE_UPPER, Xmid); + set_variable ("sigmoid_right_center", par[0]); + set_variable ("sigmoid_right_sigma", par[1]); + set_variable ("sigmoid_right_top", par[2]); + set_variable ("sigmoid_right_bottom", par[3]); + + return TRUE; + +usage: + gprint (GP_ERR, "USAGE: deimos fitprofile xprof fprof wprof Xs Xe\n"); + return FALSE; +} + +/* pars: x_o, sigma, top, bottom */ +opihi_flt rsigmoidDeimos (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) { + OHANA_UNUSED_PARAM(Npar); + + double z = (x - par[0])/par[1]; + double r = exp (z); + double s = 1 + r; + double q = 1 / s; + double f = par[2]*q + par[3]; + + // df/dp0 = df/dq * dq/dr * dr/dz * dz/dp0 + // df/dq = par[2] + // dq/ds = -s^-2 = -q/s + // ds/dr = 1 + // dr/dz = r + // dz/dp0 = -1 / par[1] + // dz/dp1 = -z / par[1] + + // df/dp0 = par[2]*(-q/s)*r*(-1/par[1]) = par[2]*q*r/(s*par[1]) + // df/dp1 = par[2]*(-q/s)*r*(-z/par[1]) = par[2]*q*r*z/(s*par[1]) + + if (dpar) { + dpar[0] = par[2]*q*r/(s*par[1]); + dpar[1] = dpar[0]*z; + dpar[2] = q; + dpar[3] = 1; + } + + return (f); +} + +/* pars: x_o, sigma, top, bottom */ +opihi_flt fsigmoidDeimos (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) { + OHANA_UNUSED_PARAM(Npar); + + double z = (x - par[0])/par[1]; + double r = exp (-z); + double s = 1 + r; + double q = 1 / s; + double f = par[2]*q + par[3]; + + // df/dp0 = df/dq * dq/dr * dr/dz * dz/dp0 + // df/dq = par[2] + // dq/ds = -s^-2 = -q/s + // ds/dr = 1 + // dr/dz = -r + // dz/dp0 = -1 / par[1] + // dz/dp1 = -z / par[1] + + // df/dp0 = -par[2]*(-q/s)*r*(-1/par[1]) = par[2]*q*r/(s*par[1]) + // df/dp1 = -par[2]*(-q/s)*r*(-z/par[1]) = par[2]*q*r*z/(s*par[1]) + + if (dpar) { + dpar[0] = -par[2]*q*r/(s*par[1]); + dpar[1] = dpar[0]*z; + dpar[2] = q; + dpar[3] = 1; + } + + return (f); +} + +# define FIT_TOLERANCE 1e-4 +# define FLT_TOLERANCE 1e-6 + +// wprofile is the true poisson error +void sigmoid_fit_irls (opihi_flt *xprofile, opihi_flt *fprofile, opihi_flt *wprofile, opihi_flt *oprofile, int NptsAll, opihi_flt *par, int Npar, SigmoidMode mode, double limit) { + + FitFunc *FUNC = (mode == MODE_LOWER) ? fsigmoidDeimos : rsigmoidDeimos; + + ALLOCATE_PTR (par_old, opihi_flt, Npar); + + // allocate these as global static in this file? + ALLOCATE_PTR (xtemp, opihi_flt, NptsAll); + ALLOCATE_PTR (ytemp, opihi_flt, NptsAll); + ALLOCATE_PTR (wtemp, opihi_flt, NptsAll); + ALLOCATE_PTR (otemp, opihi_flt, NptsAll); + + // extract the points on the right side and fit + int Npts = 0; + for (int i = 0; i < NptsAll; i++) { + if (!isfinite(wprofile[i])) continue; + if (!isfinite(fprofile[i])) continue; + int keep = (mode == MODE_UPPER) ^ (xprofile[i] < limit); + if (keep) { + xtemp[Npts] = xprofile[i]; + ytemp[Npts] = fprofile[i]; + wtemp[Npts] = (wprofile[i] == 0.0) ? 1.0 : wprofile[i]; // crude + Npts ++; + } + } + + // define initial weights as 1.0 + ALLOCATE_PTR (wt, opihi_flt, Npts); + for (int i = 0; i < Npts; i++) { + wt[i] = 1.0 / SQ(wtemp[i]); + } + + // initial fit with ordinary least-squares + // otemp is the fitted value of the function at xtemp + sigmoid_fit (xtemp, ytemp, wt, otemp, Npts, par, Npar, FUNC); + + int converged = FALSE; + for (int iterations = 0; !converged && (iterations < MaxIterations); iterations++) { + + for (int i = 0; i < Npts; i++) { + wt[i] = deimos_weight_cauchy ((ytemp[i] - otemp[i]) / wtemp[i]) / SQ(wtemp[i]); + } + + sigmoid_fit (xtemp, ytemp, wt, otemp, Npts, par, Npar, FUNC); + + // save this solution + for (int i = 0; i < Npar; i++) { + par_old[i] = par[i]; + } + + converged = TRUE; + for (int i = 0; i < Npts; i++) { + if ((fabs(par[i] - par_old[i]) > FIT_TOLERANCE * fabs(par[i])) && + (fabs(par[i] - par_old[i]) > FLT_TOLERANCE)) + converged = FALSE; + } + } + + for (int i = 0; i < NptsAll; i++) { + oprofile[i] = FUNC (xprofile[i], par, Npar, NULL); + } + + FREE (xtemp); + FREE (ytemp); + FREE (wtemp); + FREE (otemp); + + FREE (par_old); + FREE (wt); +} + +void sigmoid_fit (opihi_flt *xprofile, opihi_flt *yprofile, opihi_flt *wprofile, opihi_flt *oprofile, int Npts, opihi_flt *par, int Npar, FitFunc *FUNC) { + + double ochisq = mrqinit (xprofile, yprofile, wprofile, Npts, par, Npar, FUNC, VERBOSE); + double dchisq = ochisq + 2*Npts; + + int Niter; + for (Niter = 0; (Niter < 20) && ((dchisq > 0.1*(Npts - Npar)) || (dchisq <= 0.0)); Niter++) { + double chisq = mrqmin (xprofile, yprofile, wprofile, Npts, par, Npar, FUNC, VERBOSE); + dchisq = ochisq - chisq; + ochisq = chisq; + if (VERBOSE) gprint (GP_ERR, "dchisq: %f, Ndof: %d\n", dchisq, Npts - Npar); + } + if (VERBOSE) gprint (GP_ERR, "%d iterations\n", Niter); + + if (!QUIET) gprint (GP_LOG, "sigmoid @ %f (%f) : top %f bottom %f\n", par[0], par[1], par[2], par[3]); + + for (int i = 0; i < Npts; i++) { + oprofile[i] = FUNC (xprofile[i], par, Npar, NULL); + } + + mrqfree (Npar); +} + +double deimos_weight_cauchy (double x) { + double r = x / 2.385; + return (1.0 / (1.0 + SQ(r))); +} + Index: trunk/Ohana/src/opihi/mana/deimos_fitslit.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_fitslit.c (revision 41158) +++ trunk/Ohana/src/opihi/mana/deimos_fitslit.c (revision 41159) @@ -1,9 +1,11 @@ # include "data.h" +# define IRLS_TOLERANCE 1e-4 + +float weight_cauchy_square_flt (float x2); +static void fitflux (float *Fwind, float *Fprof, float *weight, int Nx, double *flux, double *sky); int deimos_fitslit (int argc, char **argv) { - // fitslit (window) (slit) (flux) (sky) - - // int N; + // fitslit (window) (vari) (slit) (flux) (sky) Vector *flux = NULL; @@ -11,9 +13,25 @@ Buffer *wind = NULL; + Buffer *vari = NULL; Buffer *slit = NULL; - if (argc != 5) { - gprint (GP_ERR, "USAGE: deimos fitslit (window) (slit) (flux) (sky)\n"); - gprint (GP_ERR, " inputs: window (observed 2D flux), slit (model 2D flux)\n"); + int N; + int Niter = 10; + if ((N = get_argument (argc, argv, "-irls-iter"))) { + remove_argument (N, &argc, argv); + Niter = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + Buffer *model = NULL; + if ((N = get_argument (argc, argv, "-model"))) { + remove_argument (N, &argc, argv); + if ((model = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } + + if (argc != 6) { + gprint (GP_ERR, "USAGE: deimos fitslit (window) (variance) (slit) (flux) (sky)\n"); + gprint (GP_ERR, " inputs: window (observed 2D flux), variance (on 2D flux), slit (model 2D flux)\n"); gprint (GP_ERR, " outputs: flux (best-fit 1D flux), sky (best-fit 1D background)\n"); return FALSE; @@ -22,9 +40,8 @@ // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior if ((wind = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); - if ((slit = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE); - if ((flux = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE); - if ((sky = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE); - - // XXX confirm slit and wind have same dimensions + if ((vari = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE); + if ((slit = SelectBuffer (argv[3], OLDBUFFER, TRUE)) == NULL) return (FALSE); + if ((flux = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((sky = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE); // define the output window @@ -32,42 +49,117 @@ int Ny = wind[0].matrix.Naxis[1]; + // confirm wind, vari, slit have same dimensions if (Nx != slit[0].matrix.Naxis[0]) { gprint (GP_ERR, "size mismatch in fitslit for window and slit\n"); return FALSE; } if (Ny != slit[0].matrix.Naxis[1]) { gprint (GP_ERR, "size mismatch in fitslit for window and slit\n"); return FALSE; } + if (Nx != vari[0].matrix.Naxis[0]) { gprint (GP_ERR, "size mismatch in fitslit for window and variance\n"); return FALSE; } + if (Ny != vari[0].matrix.Naxis[1]) { gprint (GP_ERR, "size mismatch in fitslit for window and variance\n"); return FALSE; } ResetVector (flux, OPIHI_FLT, Ny); ResetVector (sky, OPIHI_FLT, Ny); + if (model) { + if (!CreateBuffer (model, Nx, Ny, -32, 1.0, 0.0)) { gprint (GP_ERR, "error generating output model buffer\n"); return FALSE; } + } + float *Fwind = (float *) wind[0].matrix.buffer; + float *Fvari = (float *) vari[0].matrix.buffer; float *Fprof = (float *) slit[0].matrix.buffer; + float *FmodOut = model ? (float *) model[0].matrix.buffer : NULL; + + ALLOCATE_PTR (weight, float, Nx); + ALLOCATE_PTR (rawwgt, float, Nx); // loop over the rows for (int iy = 0; iy < Ny; iy++) { - // calculate elements of the chi-square for this row: + double F = 0.0, S = 0.0; - double R = 0.0, P = 0.0, P2 = 0.0, f = 0.0, fp = 0.0; - double wt = 1.0; - + // set weight based on the variance for (int ix = 0; ix < Nx; ix++) { - - float fxy = Fwind[ix + iy*Nx]; - float pxy = Fprof[ix + iy*Nx]; - - R += wt; - P += pxy*wt; - P2 += pxy*pxy*wt; - f += fxy*wt; - fp += fxy*pxy*wt; + float vari = Fvari[ix + iy*Nx]; + weight[ix] = !isfinite(vari) || (fabs(vari) < 1e-6) ? 1.0 : 1.0 / vari; + rawwgt[ix] = weight[ix]; } - double det = R*P2 - P*P; - double S = (P2*f - P*fp) / det; - double F = (R*fp - P*f) / det; + // calculate elements of the chi-square for this row: + fitflux (&Fwind[iy*Nx], &Fprof[iy*Nx], weight, Nx, &F, &S); + + // do an IRLS loop here to reject outliers + int converged = FALSE; + for (int iter = 0; (iter < Niter) && !converged; iter++) { + + double Flast = F; + double Slast = S; + + // calculate weight modification based on distances (squared). + // use modifier to calculate new weighted mean + for (int ix = 0; ix < Nx; ix++) { + + // F & S are model parameters; compare the model and observed values for each pixel + float Fmodel = S + F*Fprof[ix + iy*Nx]; + float dV = (Fwind[ix + iy*Nx] - Fmodel); + float d2 = SQ(dV) * rawwgt[ix]; + + float Mod = weight_cauchy_square_flt (d2); + weight[ix] = Mod * rawwgt[ix]; + } + fitflux (&Fwind[iy*Nx], &Fprof[iy*Nx], weight, Nx, &F, &S); + + float dF = fabs(F - Flast); + float dS = fabs(S - Slast); + + if ((dF < F * IRLS_TOLERANCE) && (dS < S * IRLS_TOLERANCE)) converged = TRUE; + } + + if (FmodOut) { + for (int ix = 0; ix < Nx; ix++) { + // F & S are model parameters; compare the model and observed values for each pixel + float Fmodel = S + F*Fprof[ix + iy*Nx]; + float dV = (Fwind[ix + iy*Nx] - Fmodel); + float d2 = SQ(dV) * rawwgt[ix]; + float Mod = weight_cauchy_square_flt (d2); + if (Mod < 0.1) { + FmodOut[ix + iy*Nx] = NAN; + } else { + FmodOut[ix + iy*Nx] = Fmodel; + } + } + } flux->elements.Flt[iy] = F; sky->elements.Flt[iy] = S; } - return TRUE; } +// weight is 1 / variance +// Fwind, Fprof are pointers to the start of a single row +// Fwind is the observed flux in the window +// Fprof is the slit profile +void fitflux (float *Fwind, float *Fprof, float *weight, int Nx, double *flux, double *sky) { + + // calculate elements of the chi-square fit for this row: + + double R = 0.0, P = 0.0, P2 = 0.0, f = 0.0, fp = 0.0; + + for (int ix = 0; ix < Nx; ix++) { + + double fxy = Fwind[ix]; + double pxy = Fprof[ix]; + + // skip NAN / Inf pixels in either window or slit + if (!isfinite(fxy)) continue; + if (!isfinite(pxy)) continue; + + double wt = (weight == NULL) ? 1.0 : weight[ix]; + R += wt; + P += pxy*wt; + P2 += pxy*pxy*wt; + f += fxy*wt; + fp += fxy*pxy*wt; + } + + double det = R*P2 - P*P; + *sky = (P2*f - P*fp) / det; + *flux = (R*fp - P*f) / det; +} Index: trunk/Ohana/src/opihi/mana/deimos_getobj.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_getobj.c (revision 41158) +++ trunk/Ohana/src/opihi/mana/deimos_getobj.c (revision 41159) @@ -3,29 +3,17 @@ // temporary storage used in local functions: -float *tmpv = NULL; -float *tmpx = NULL; - -// pass in a pointer to the row (&buffV[iy*Nx]) -float get_segment_median (float *buffV, float Xs, float Xe, float *Xref) { - - // get the pixels in the left portion below profile: - // XXX: worry about the effect of a fractional offset - for (int ix = Xs + 0.5; ix < Xe + 0.5; ix++) { - int i = ix - Xs; - tmpv[i] = buffV[ix]; - tmpx[i] = ix; - } - int Nv = Xe - Xs; - - // generate a statistic from these pixels: mean, median, etc? - fsort (tmpv, Nv); - int Nv2 = Nv/2; // Nv == 1, Nv2 = 0, Nv == 2, Nv2 = 1, Nv == 3, Nv2 = 1, Nv == 4, Nv == 2 - float V = (Nv % 2) ? tmpv[Nv2] : 0.5*(tmpv[Nv2-1] + tmpv[Nv2]); - - // calculate average coordinate (not flux-weighted...) - // XXX : if we allow (and exclude) NAN-valued pixels, we can calculate this more explicitly - *Xref = 0.5*(Xe + Xs); - - return V; +static float *tmpv = NULL; +static float *tmpx = NULL; + +float get_median_flt (float *val, int Nval) { + + if (Nval <= 0) return NAN; + + int Nval2 = Nval / 2; + + fsort (val, Nval); + float value = (Nval % 2) ? val[Nval2] : 0.5*(val[Nval2-1] + val[Nval2]); + + return value; } @@ -37,7 +25,12 @@ // not used : LSF -- not needed (measurement per disp pixel) + // This version uses the median of the pixels in the image which are + // in regions where the slit profile is < 15% (as opposed to the + // getobj version which attempts to fit both background and sky regions + // input parameters: // * buffer : 2D image of slit (full frame or cutout?) - // * trace : spline fit of slit central x pos vs y-coord + // * slit_trace : spline fit of slit central x pos vs y-coord + // * psf_trace : spline fit of psf central x pos vs y-coord // * profile : slit window profile (vector) // * PSF : point-spread function vector (flux normalized, x-dir) @@ -51,5 +44,7 @@ int N; - Spline *trace = NULL; + Spline *slit_trace_red = NULL; + Spline *slit_trace_blu = NULL; + Spline *psf_trace = NULL; Vector *profile = NULL; @@ -64,14 +59,28 @@ float stilt = 0.0; // angle of the slit - int SaveSegments = FALSE; - if ((N = get_argument (argc, argv, "-save-segments"))) { - remove_argument (N, &argc, argv); - SaveSegments = TRUE; - } + // for a red vs blu spline, we need to specify the split point + // XXX this is REALLY ad-hoc for Deimos. not sure how to make + // this more generic (need to define the ranges somewhere) + int redlimit = 4096; + if ((N = get_argument (argc, argv, "-redlimit"))) { + remove_argument (N, &argc, argv); + redlimit = atoi (argv[N]); + remove_argument (N, &argc, argv); + } // Input parameters: - if ((N = get_argument (argc, argv, "-trace"))) { - remove_argument (N, &argc, argv); - if ((trace = FindSpline (argv[N])) == NULL) return (FALSE); + if ((N = get_argument (argc, argv, "-slit-trace-red"))) { + remove_argument (N, &argc, argv); + if ((slit_trace_red = FindSpline (argv[N])) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-slit-trace-blu"))) { + remove_argument (N, &argc, argv); + if ((slit_trace_blu = FindSpline (argv[N])) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-psf-trace"))) { + remove_argument (N, &argc, argv); + if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE); remove_argument (N, &argc, argv); } else { goto usage; } @@ -137,47 +146,4 @@ ResetVector (backgnd, OPIHI_FLT, Ny); - // I need to generate a sky-left and sky-right vector in a region on the slit portion of - // the profile, but off the object. Then each of these vectors can be shifted up and down - // based on stilt. Use these sky vectors to predict the sky under the PSF - - // *** scan the profile to determine the background and sky cut points *** - - // values below in the profile coordinate system need to be shifted by (Xref - Nprofile/2) - // to match the image coordinate system - - int B_s_l = 0, B_e_l = -1; // left background window - int B_s_r = -1, B_e_r = Nprofile - 1; // right background window - int S_left = -1, S_right = -1; // top portion of profile - for (int i = 0; i < Nprofile; i++) { - if ((B_e_l < 0) && (profile->elements.Flt[i] > 0.10)) B_e_l = i - 1; - if ((B_e_l >= 0) && (S_left < 0) && (profile->elements.Flt[i] > 0.95)) S_left = i; - if ((S_left >= 0) && (S_right < 0) && (profile->elements.Flt[i] < 0.95)) S_right = i - 1; - if ((S_right >= 0) && (profile->elements.Flt[i] < 0.10)) { B_s_r = i; break; } - } - if (B_e_l == B_s_l) B_e_l ++; - if (B_e_r == B_s_r) B_s_r --; - - // S_left to S_right defines are window at the top, choose 0.05 - 0.15, 0.85 - 0.95 - int N_top = S_right - S_left + 1; - - // if N_top < 20, these jumps might be too small - // if N_top < 20, these jumps might be too small - int S_s_l = S_left + 0.05*N_top; - int S_e_l = S_left + 0.15*N_top; - if (S_e_l == S_s_l) S_e_l ++; - int S_s_r = S_left + 0.85*N_top; - int S_e_r = S_left + 0.95*N_top; - if (S_e_r == S_s_r) S_s_r --; - - gprint (GP_ERR, "windows: %.1f - %.1f, %.1f - %.1f, %.1f - %.1f, %.1f - %.1f\n", - B_s_l + Xref - Xprofile, B_e_l + Xref - Xprofile, - S_s_l + Xref - Xprofile, S_e_l + Xref - Xprofile, - S_s_r + Xref - Xprofile, S_e_r + Xref - Xprofile, - B_s_r + Xref - Xprofile, B_e_r + Xref - Xprofile); - - // OPEN QUESTIONS: - // NOTE: if the range of pixels in the windows is too large, then the tilt will wash out the signal. - // this pass is a guess, so perhaps that does not matter, but there is clearly a trade-off here. - // generate two buffers to store the flux and coordinate for the background, sky, object pixels: ALLOCATE (tmpv, float, Nx); @@ -186,75 +152,9 @@ float *buffV = (float *) buffer[0].matrix.buffer; - // use these to store the flux vectors - ALLOCATE_PTR (bckF_lf, float, Ny); - ALLOCATE_PTR (bckF_rt, float, Ny); - ALLOCATE_PTR (skyF_lf, float, Ny); - ALLOCATE_PTR (skyF_rt, float, Ny); - - // use these to store the x-coord vectors - ALLOCATE_PTR (bckX_lf, float, Ny); - ALLOCATE_PTR (bckX_rt, float, Ny); - ALLOCATE_PTR (skyX_lf, float, Ny); - ALLOCATE_PTR (skyX_rt, float, Ny); - float dxMax = 0; // maximum distance between sky regions - Vector *skyLsV = NULL, *skyLeV = NULL, *skyLfV = NULL, *skyRsV = NULL, *skyReV = NULL, *skyRfV = NULL; - Vector *bckLsV = NULL, *bckLeV = NULL, *bckLfV = NULL, *bckRsV = NULL, *bckReV = NULL, *bckRfV = NULL; - if (SaveSegments) { - skyLsV = SelectVector ("skyLs", ANYVECTOR, TRUE); ResetVector (skyLsV, OPIHI_FLT, Ny); - skyLeV = SelectVector ("skyLe", ANYVECTOR, TRUE); ResetVector (skyLeV, OPIHI_FLT, Ny); - skyLfV = SelectVector ("skyLf", ANYVECTOR, TRUE); ResetVector (skyLfV, OPIHI_FLT, Ny); - skyRsV = SelectVector ("skyRs", ANYVECTOR, TRUE); ResetVector (skyRsV, OPIHI_FLT, Ny); - skyReV = SelectVector ("skyRe", ANYVECTOR, TRUE); ResetVector (skyReV, OPIHI_FLT, Ny); - skyRfV = SelectVector ("skyRf", ANYVECTOR, TRUE); ResetVector (skyRfV, OPIHI_FLT, Ny); - bckLsV = SelectVector ("bckLs", ANYVECTOR, TRUE); ResetVector (bckLsV, OPIHI_FLT, Ny); - bckLeV = SelectVector ("bckLe", ANYVECTOR, TRUE); ResetVector (bckLeV, OPIHI_FLT, Ny); - bckLfV = SelectVector ("bckLf", ANYVECTOR, TRUE); ResetVector (bckLfV, OPIHI_FLT, Ny); - bckRsV = SelectVector ("bckRs", ANYVECTOR, TRUE); ResetVector (bckRsV, OPIHI_FLT, Ny); - bckReV = SelectVector ("bckRe", ANYVECTOR, TRUE); ResetVector (bckReV, OPIHI_FLT, Ny); - bckRfV = SelectVector ("bckRf", ANYVECTOR, TRUE); ResetVector (bckRfV, OPIHI_FLT, Ny); - } - - // loop over the y positions - for (int iy = 0; iy < Ny; iy++) { - - // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center - float dXtrace = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy); - - // Xref is the reference coordinate of the profile on the image - // Xprofile is the center of the profile. - // dXtrace follows the x-dir variations of the slit profile - float Xtrace = Xref - Xprofile + dXtrace; - - bckF_lf[iy] = get_segment_median(&buffV[iy*Nx], B_s_l + Xtrace, B_e_l + Xtrace, &bckX_lf[iy]); - bckF_rt[iy] = get_segment_median(&buffV[iy*Nx], B_s_r + Xtrace, B_e_r + Xtrace, &bckX_rt[iy]); - skyF_lf[iy] = get_segment_median(&buffV[iy*Nx], S_s_l + Xtrace, S_e_l + Xtrace, &skyX_lf[iy]); - skyF_rt[iy] = get_segment_median(&buffV[iy*Nx], S_s_r + Xtrace, S_e_r + Xtrace, &skyX_rt[iy]); - dxMax = MAX (dxMax, skyX_rt[iy] - skyX_lf[iy]); // this is a constant value unless we ignore some of the pixels - // assert skyX_rt > skyX_lf ? - - if (SaveSegments) { - skyLsV->elements.Flt[iy] = S_s_l + Xtrace; - skyLeV->elements.Flt[iy] = S_e_l + Xtrace; - skyLfV->elements.Flt[iy] = skyF_lf[iy]; - skyRsV->elements.Flt[iy] = S_s_r + Xtrace; - skyReV->elements.Flt[iy] = S_e_r + Xtrace; - skyRfV->elements.Flt[iy] = skyF_rt[iy]; - bckLsV->elements.Flt[iy] = B_s_l + Xtrace; - bckLeV->elements.Flt[iy] = B_e_l + Xtrace; - bckLfV->elements.Flt[iy] = bckF_lf[iy]; - bckRsV->elements.Flt[iy] = B_s_r + Xtrace; - bckReV->elements.Flt[iy] = B_e_r + Xtrace; - bckRfV->elements.Flt[iy] = bckF_rt[iy]; - } - } - // extract sky-subtracted values for the psf fit, using the two sky vectors to predict the sky - float dyMax = dxMax*sin(stilt*RAD_DEG); // dxMax is based on the profile width, calculated above - - // find nominal dy_lf, dy_rt offsets for the sky vectors. - int dy_lf = -dyMax / 2; - int dy_rt = +dyMax / 2; + float sin_stilt = sin(stilt*RAD_DEG); + float dyMax = dxMax*sin_stilt; // dxMax is based on the profile width, calculated above // we are fitting this PSF to the data: @@ -273,32 +173,99 @@ } - - // Ys, Ye based on the dyMax value above (skip the first few lines to avoid running out of bounds on sky) for (int iy = fabs(dyMax); iy < Ny - fabs(dyMax); iy++) { - // X_lf & X_rt are calculated in image x-coords - float S_lf = skyF_lf[iy + dy_lf]; - float S_rt = skyF_rt[iy + dy_rt]; - float X_lf = skyX_lf[iy + dy_lf]; - float X_rt = skyX_rt[iy + dy_rt]; - // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center - float dXtrace = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy); + Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu; + float dXtraceSLIT = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy); + float dXtracePSF = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy); // Xref is the reference coordinate of the profile on the image // Xpsf is the center of the PSF. // dXtrace follows the x-dir variations of the slit profile (or PSF?) - float Xtrace = Xref - Xpsf + dXtrace; - - // include the PSF trace and restrict the fit below to portions where the profile is > 0.95 + + // XtracePSF is the offset between the buffer ix and the psf elements at this location + float XtracePSF = Xref - Xpsf + dXtraceSLIT + dXtracePSF; + + // XtraceSLIT is the offset between the buffer ix and the slit profile elements at this location + float XtraceSLIT = Xref - Xprofile + dXtraceSLIT; + + // measure the background outside the slit + + // select pixels from the image which are in the background regions (profile < 0.15) + int Nbck = 0; + int Nsky = 0; + for (int ix = 0; ix < Nx; ix ++) { + int islit = ix - XtraceSLIT - 0.5; // slit window coordinate + int ipsf = ix - XtracePSF - 0.5; // psf model coordinate + + if (islit < 0) continue; + if (islit >= Nprofile) continue; + + // coordinate of the PSF center: + int dxpsf = ipsf - Xpsf; + int dypsf = dxpsf * sin_stilt; + + int iY = iy + dypsf; + + if (iy == 6000) { + fprintf (stderr, "bck : %d %d %d : %f : %f\n", ix, islit, ipsf, buffV[ix + iY*Nx], profile->elements.Flt[islit]); + } + + // pixels in the off-slit background: + if (profile->elements.Flt[islit] < 0.15) { + tmpv[Nbck] = buffV[ix + iY*Nx]; + if (!isfinite(tmpv[Nbck])) continue; + Nbck ++; + continue; + } + + // the psf window can be smaller than the slit window + if (ipsf < 0) continue; + if (ipsf >= Npsf) continue; + + // pixels in the on-slit sky: + if ((profile->elements.Flt[islit] > 0.95) && (psfV[ipsf] < 0.15)) { + tmpx[Nsky] = buffV[ix + iY*Nx]; + if (!isfinite(tmpx[Nsky])) continue; + Nsky ++; + } + } + backgnd->elements.Flt[iy] = get_median_flt (tmpv, Nbck); + float So = get_median_flt (tmpx, Nsky); + sky->elements.Flt[iy] = So; // now I need to find the limits of the PSF in the buffer at this iy location - for (int ipsf = 0; ipsf < Npsf; ipsf ++) { - - int ix = ipsf + Xtrace + 0.5; // image coordinate (nearest pixel center) - float So = S_lf + (ix - X_lf)*(S_rt - S_lf)/(X_rt - X_lf); // interpolated sky under PSF peak - - tmpv[ipsf] = buffV[ix + iy*Nx] - So; - tmpx[ipsf] = ipsf; + int Ntmp = 0; + for (int ix = 0; ix < Nx; ix ++) { + int islit = ix - XtraceSLIT - 0.5; // slit window coordinate + int ipsf = ix - XtracePSF - 0.5; // psf model coordinate + + if (islit < 0) continue; + if (islit >= Nprofile) continue; + + // the psf window can be smaller than the slit window + if (ipsf < 0) continue; + if (ipsf >= Npsf) continue; + + // coordinate of the PSF center: + int dxpsf = ipsf - Xpsf; + int dypsf = dxpsf * sin_stilt; + + int iY = iy + dypsf; + + // filter on the PSF value? + if (profile->elements.Flt[islit] < 0.95) continue; + + if (iy == 6000) { + fprintf (stderr, "%d %d %d : %f %f : %f : %f\n", ix, islit, ipsf, buffV[ix + iY*Nx] - So, psfV[ipsf], (buffV[ix + iY*Nx] - So)*psfV[ipsf], profile->elements.Flt[islit]); + } + + tmpv[Ntmp] = buffV[ix + iY*Nx] - So; + tmpx[Ntmp] = psfV[ipsf]; + + if (!isfinite(tmpv[Ntmp])) continue; + if (!isfinite(tmpx[Ntmp])) continue; + + Ntmp ++; } @@ -306,31 +273,16 @@ // now fit the PSF to tmpv by calculating a normalization float S1 = 0.0, S2 = 0.0; - for (int ipsf = 0; ipsf < Npsf; ipsf++) { + for (int ipsf = 0; ipsf < Ntmp; ipsf++) { if (iy == 6000) { - fprintf (stderr, "%d : %f %f : %f\n", ipsf, tmpv[ipsf], psfV[ipsf], tmpv[ipsf]*psfV[ipsf]); - } - S1 += tmpv[ipsf]*psfV[ipsf]; // if we have an errorbar, include / SQ(sigv[i]) - S2 += psfV[ipsf]*psfV[ipsf]; // if we have an errorbar, include / SQ(sigv[i]) + fprintf (stderr, "%d : %f %f : %f\n", ipsf, tmpv[ipsf], tmpx[ipsf], tmpv[ipsf]*tmpx[ipsf]); + } + S1 += tmpv[ipsf]*tmpx[ipsf]; // if we have an errorbar, include / SQ(sigv[i]) + S2 += tmpx[ipsf]*tmpx[ipsf]; // if we have an errorbar, include / SQ(sigv[i]) } - - int Xo = Xref - dXtrace; // center of PSF in image x-coordinates - float So = S_lf + (Xo - X_lf)*(S_rt - S_lf)/(X_rt - X_lf); // interpolated sky under PSF peak object->elements.Flt[iy] = S1 / S2; - sky->elements.Flt[iy] = So; - backgnd->elements.Flt[iy] = 0.5*(bckF_lf[iy] + bckF_rt[iy]); } free (tmpv); free (tmpx); - - free (bckF_lf); - free (bckF_rt); - free (skyF_lf); - free (skyF_rt); - - free (bckX_lf); - free (bckX_rt); - free (skyX_lf); - free (skyX_rt); return TRUE; Index: trunk/Ohana/src/opihi/mana/deimos_getobj.v0.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_getobj.v0.c (revision 41159) +++ trunk/Ohana/src/opihi/mana/deimos_getobj.v0.c (revision 41159) @@ -0,0 +1,383 @@ +# include "data.h" + +// temporary storage used in local functions: + +static float *tmpv = NULL; +static float *tmpx = NULL; + +// pass in a pointer to the row (&buffV[iy*Nx]) +float get_segment_median (float *buffV, float Xs, float Xe, float *Xref) { + + // get the pixels in the left portion below profile: + // XXX: worry about the effect of a fractional offset + for (int ix = Xs + 0.5; ix < Xe + 0.5; ix++) { + int i = ix - Xs; + tmpv[i] = buffV[ix]; + tmpx[i] = ix; + } + int Nv = Xe - Xs; + + // generate a statistic from these pixels: mean, median, etc? + fsort (tmpv, Nv); + int Nv2 = Nv/2; // Nv == 1, Nv2 = 0, Nv == 2, Nv2 = 1, Nv == 3, Nv2 = 1, Nv == 4, Nv == 2 + float V = (Nv % 2) ? tmpv[Nv2] : 0.5*(tmpv[Nv2-1] + tmpv[Nv2]); + + // calculate average coordinate (not flux-weighted...) + // XXX : if we allow (and exclude) NAN-valued pixels, we can calculate this more explicitly + *Xref = 0.5*(Xe + Xs); + + return V; +} + +int deimos_getobj (int argc, char **argv) { + + // determine rough object model parameters based on quick stats + // input : buffer PSF trace profile stilt + // output : object sky background + // not used : LSF -- not needed (measurement per disp pixel) + + // input parameters: + // * buffer : 2D image of slit (full frame or cutout?) + // * slit_trace : spline fit of slit central x pos vs y-coord + // * psf_trace : spline fit of psf central x pos vs y-coord + // * profile : slit window profile (vector) + // * PSF : point-spread function vector (flux normalized, x-dir) + // * stilt : slit tilt response : 2D kernel? + + // output parameters: + // * object : vector of object flux vs y-coord + // * sky : vector of local sky signal vs y-coord + // * background : vector of extra-slit background flux vs y-coord + + int N; + + Spline *slit_trace = NULL; + Spline *psf_trace = NULL; + Vector *profile = NULL; + + Vector *object = NULL; + Vector *sky = NULL; + Vector *backgnd = NULL; + + Vector *psf = NULL; + + Buffer *buffer = NULL; + + float stilt = 0.0; // angle of the slit + + int SaveSegments = FALSE; + if ((N = get_argument (argc, argv, "-save-segments"))) { + remove_argument (N, &argc, argv); + SaveSegments = TRUE; + } + + // Input parameters: + if ((N = get_argument (argc, argv, "-slit-trace"))) { + remove_argument (N, &argc, argv); + if ((slit_trace = FindSpline (argv[N])) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-psf-trace"))) { + remove_argument (N, &argc, argv); + if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-profile"))) { + remove_argument (N, &argc, argv); + if ((profile = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-psf"))) { + remove_argument (N, &argc, argv); + if ((psf = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-stilt"))) { + remove_argument (N, &argc, argv); + stilt = atof (argv[N]); + remove_argument (N, &argc, argv); + } else { goto usage; } + + // Output parameters: + if ((N = get_argument (argc, argv, "-object"))) { + remove_argument (N, &argc, argv); + if ((object = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-sky"))) { + remove_argument (N, &argc, argv); + if ((sky = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-backgnd"))) { + remove_argument (N, &argc, argv); + if ((backgnd = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + + if (argc != 3) goto usage; + + // buffer is the full image, Xref is reference coordinate of the profile in the buffer + if ((buffer = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + float Xref = atof(argv[2]); + + // profile and PSF are defined with central reference pixel mapped to Xref (+ trace) + int Nprofile = profile->Nelements; + if (Nprofile % 2 == 0) { + gprint (GP_ERR, "slit profile vector must have an odd number of pixels\n"); + return FALSE; + } + int Xprofile = Nprofile / 2; + + int Npsf = psf->Nelements; + if (Npsf % 2 == 0) { + gprint (GP_ERR, "PSF vector must have an odd number of pixels\n"); + return FALSE; + } + int Xpsf = Npsf / 2; + + // we will generate output vectors (object, sky, background) of length Ny + int Nx = buffer[0].matrix.Naxis[0]; + int Ny = buffer[0].matrix.Naxis[1]; + ResetVector (object, OPIHI_FLT, Ny); + ResetVector (sky, OPIHI_FLT, Ny); + ResetVector (backgnd, OPIHI_FLT, Ny); + + // I need to generate a sky-left and sky-right vector in a region on the slit portion of + // the profile, but off the object. Then each of these vectors can be shifted up and down + // based on stilt. Use these sky vectors to predict the sky under the PSF + + // *** scan the profile to determine the background and sky cut points *** + + // values below in the profile coordinate system need to be shifted by (Xref - Nprofile/2) + // to match the image coordinate system + + int B_s_l = 0, B_e_l = -1; // left background window + int B_s_r = -1, B_e_r = Nprofile - 1; // right background window + int S_left = -1, S_right = -1; // top portion of profile + for (int i = 0; i < Nprofile; i++) { + if ((B_e_l < 0) && (profile->elements.Flt[i] > 0.10)) B_e_l = i - 1; + if ((B_e_l >= 0) && (S_left < 0) && (profile->elements.Flt[i] > 0.95)) S_left = i; + if ((S_left >= 0) && (S_right < 0) && (profile->elements.Flt[i] < 0.95)) S_right = i - 1; + if ((S_right >= 0) && (profile->elements.Flt[i] < 0.10)) { B_s_r = i; break; } + } + if (B_e_l == B_s_l) B_e_l ++; + if (B_e_r == B_s_r) B_s_r --; + + // S_left to S_right defines are window at the top, choose 0.05 - 0.15, 0.85 - 0.95 + int N_top = S_right - S_left + 1; + + // if N_top < 20, these jumps might be too small + // if N_top < 20, these jumps might be too small + int S_s_l = S_left + 0.05*N_top; + int S_e_l = S_left + 0.15*N_top; + if (S_e_l == S_s_l) S_e_l ++; + int S_s_r = S_left + 0.85*N_top; + int S_e_r = S_left + 0.95*N_top; + if (S_e_r == S_s_r) S_s_r --; + + gprint (GP_ERR, "windows: %.1f - %.1f, %.1f - %.1f, %.1f - %.1f, %.1f - %.1f\n", + B_s_l + Xref - Xprofile, B_e_l + Xref - Xprofile, + S_s_l + Xref - Xprofile, S_e_l + Xref - Xprofile, + S_s_r + Xref - Xprofile, S_e_r + Xref - Xprofile, + B_s_r + Xref - Xprofile, B_e_r + Xref - Xprofile); + + // OPEN QUESTIONS: + // NOTE: if the range of pixels in the windows is too large, then the tilt will wash out the signal. + // this pass is a guess, so perhaps that does not matter, but there is clearly a trade-off here. + + // generate two buffers to store the flux and coordinate for the background, sky, object pixels: + ALLOCATE (tmpv, float, Nx); + ALLOCATE (tmpx, float, Nx); + + float *buffV = (float *) buffer[0].matrix.buffer; + + // use these to store the flux vectors + ALLOCATE_PTR (bckF_lf, float, Ny); + ALLOCATE_PTR (bckF_rt, float, Ny); + ALLOCATE_PTR (skyF_lf, float, Ny); + ALLOCATE_PTR (skyF_rt, float, Ny); + + // use these to store the x-coord vectors + ALLOCATE_PTR (bckX_lf, float, Ny); + ALLOCATE_PTR (bckX_rt, float, Ny); + ALLOCATE_PTR (skyX_lf, float, Ny); + ALLOCATE_PTR (skyX_rt, float, Ny); + + float dxMax = 0; // maximum distance between sky regions + + Vector *skyLsV = NULL, *skyLeV = NULL, *skyLfV = NULL, *skyRsV = NULL, *skyReV = NULL, *skyRfV = NULL; + Vector *bckLsV = NULL, *bckLeV = NULL, *bckLfV = NULL, *bckRsV = NULL, *bckReV = NULL, *bckRfV = NULL; + if (SaveSegments) { + skyLsV = SelectVector ("skyLs", ANYVECTOR, TRUE); ResetVector (skyLsV, OPIHI_FLT, Ny); + skyLeV = SelectVector ("skyLe", ANYVECTOR, TRUE); ResetVector (skyLeV, OPIHI_FLT, Ny); + skyLfV = SelectVector ("skyLf", ANYVECTOR, TRUE); ResetVector (skyLfV, OPIHI_FLT, Ny); + skyRsV = SelectVector ("skyRs", ANYVECTOR, TRUE); ResetVector (skyRsV, OPIHI_FLT, Ny); + skyReV = SelectVector ("skyRe", ANYVECTOR, TRUE); ResetVector (skyReV, OPIHI_FLT, Ny); + skyRfV = SelectVector ("skyRf", ANYVECTOR, TRUE); ResetVector (skyRfV, OPIHI_FLT, Ny); + bckLsV = SelectVector ("bckLs", ANYVECTOR, TRUE); ResetVector (bckLsV, OPIHI_FLT, Ny); + bckLeV = SelectVector ("bckLe", ANYVECTOR, TRUE); ResetVector (bckLeV, OPIHI_FLT, Ny); + bckLfV = SelectVector ("bckLf", ANYVECTOR, TRUE); ResetVector (bckLfV, OPIHI_FLT, Ny); + bckRsV = SelectVector ("bckRs", ANYVECTOR, TRUE); ResetVector (bckRsV, OPIHI_FLT, Ny); + bckReV = SelectVector ("bckRe", ANYVECTOR, TRUE); ResetVector (bckReV, OPIHI_FLT, Ny); + bckRfV = SelectVector ("bckRf", ANYVECTOR, TRUE); ResetVector (bckRfV, OPIHI_FLT, Ny); + } + + // loop over the y positions + // examine the slit window profile to measure background and slit sky + for (int iy = 0; iy < Ny; iy++) { + + // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center + float dXtrace = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy); + + // Xref is the reference coordinate of the profile on the image + // Xprofile is the center of the profile. + // dXtrace follows the x-dir variations of the slit profile + float Xtrace = Xref - Xprofile + dXtrace; + + bckF_lf[iy] = get_segment_median(&buffV[iy*Nx], B_s_l + Xtrace, B_e_l + Xtrace, &bckX_lf[iy]); + bckF_rt[iy] = get_segment_median(&buffV[iy*Nx], B_s_r + Xtrace, B_e_r + Xtrace, &bckX_rt[iy]); + skyF_lf[iy] = get_segment_median(&buffV[iy*Nx], S_s_l + Xtrace, S_e_l + Xtrace, &skyX_lf[iy]); + skyF_rt[iy] = get_segment_median(&buffV[iy*Nx], S_s_r + Xtrace, S_e_r + Xtrace, &skyX_rt[iy]); + dxMax = MAX (dxMax, skyX_rt[iy] - skyX_lf[iy]); // this is a constant value unless we ignore some of the pixels + // assert skyX_rt > skyX_lf ? + + if (SaveSegments) { + skyLsV->elements.Flt[iy] = S_s_l + Xtrace; + skyLeV->elements.Flt[iy] = S_e_l + Xtrace; + skyLfV->elements.Flt[iy] = skyF_lf[iy]; + skyRsV->elements.Flt[iy] = S_s_r + Xtrace; + skyReV->elements.Flt[iy] = S_e_r + Xtrace; + skyRfV->elements.Flt[iy] = skyF_rt[iy]; + bckLsV->elements.Flt[iy] = B_s_l + Xtrace; + bckLeV->elements.Flt[iy] = B_e_l + Xtrace; + bckLfV->elements.Flt[iy] = bckF_lf[iy]; + bckRsV->elements.Flt[iy] = B_s_r + Xtrace; + bckReV->elements.Flt[iy] = B_e_r + Xtrace; + bckRfV->elements.Flt[iy] = bckF_rt[iy]; + } + } + + // extract sky-subtracted values for the psf fit, using the two sky vectors to predict the sky + float dyMax = dxMax*sin(stilt*RAD_DEG); // dxMax is based on the profile width, calculated above + + // find nominal dy_lf, dy_rt offsets for the sky vectors. + int dy_lf = -dyMax / 2; + int dy_rt = +dyMax / 2; + + // we are fitting this PSF to the data: + opihi_flt *psfV = psf->elements.Flt; + + // init the regions which are outside the valid region + for (int iy = 0; iy < fabs(dyMax); iy++) { + object->elements.Flt[iy] = 0; + sky->elements.Flt[iy] = 0; + backgnd->elements.Flt[iy] = 0; + } + for (int iy = Ny - fabs(dyMax); iy < Ny; iy++) { + object->elements.Flt[iy] = 0; + sky->elements.Flt[iy] = 0; + backgnd->elements.Flt[iy] = 0; + } + + + // Ys, Ye based on the dyMax value above (skip the first few lines to avoid running out of bounds on sky) + for (int iy = fabs(dyMax); iy < Ny - fabs(dyMax); iy++) { + + // X_lf & X_rt are calculated in image x-coords + float S_lf = skyF_lf[iy + dy_lf]; + float S_rt = skyF_rt[iy + dy_rt]; + float X_lf = skyX_lf[iy + dy_lf]; + float X_rt = skyX_rt[iy + dy_rt]; + + // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center + float dXtraceSLIT = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy); + float dXtracePSF = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy); + + // Xref is the reference coordinate of the profile on the image + // Xpsf is the center of the PSF. + // dXtrace follows the x-dir variations of the slit profile (or PSF?) + float XtracePSF = Xref - Xpsf + dXtraceSLIT + dXtracePSF; + float XtraceSLIT = Xref - Xprofile + dXtraceSLIT; + + // include the PSF trace and restrict the fit below to portions where the profile is > 0.95 + + // now I need to find the limits of the PSF in the buffer at this iy location + int Ntmp = 0; + for (int ipsf = 0; ipsf < Npsf; ipsf ++) { + + int ix = ipsf + XtracePSF + 0.5; // image coordinate (nearest pixel center) + int islit = ix - XtraceSLIT - 0.5; // image coordinate (nearest pixel center) + + float So = S_lf + (ix - X_lf)*(S_rt - S_lf)/(X_rt - X_lf); // interpolated sky under PSF peak + + if (islit < 0) continue; + if (islit >= profile->Nelements) continue; + + if (iy == 6000) { + fprintf (stderr, "%d : %f %f : %f : %f\n", ipsf, buffV[ix + iy*Nx] - So, psfV[ipsf], (buffV[ix + iy*Nx] - So)*psfV[ipsf], profile->elements.Flt[islit]); + } + + if (profile->elements.Flt[islit] < 0.95) continue; + + tmpv[Ntmp] = buffV[ix + iy*Nx] - So; + tmpx[Ntmp] = psfV[ipsf]; + Ntmp ++; + } + + // I now have a pair of vectors (tmpx, tmpv) which are aligned to the PSF at this location + // now fit the PSF to tmpv by calculating a normalization + float S1 = 0.0, S2 = 0.0; + for (int ipsf = 0; ipsf < Ntmp; ipsf++) { + if (iy == 6000) { + fprintf (stderr, "%d : %f %f : %f\n", ipsf, tmpv[ipsf], tmpx[ipsf], tmpv[ipsf]*tmpx[ipsf]); + } + S1 += tmpv[ipsf]*tmpx[ipsf]; // if we have an errorbar, include / SQ(sigv[i]) + S2 += tmpx[ipsf]*tmpx[ipsf]; // if we have an errorbar, include / SQ(sigv[i]) + } + + int Xo = Xref - dXtraceSLIT - dXtracePSF; // center of PSF in image x-coordinates + float So = S_lf + (Xo - X_lf)*(S_rt - S_lf)/(X_rt - X_lf); // interpolated sky under PSF peak + object->elements.Flt[iy] = S1 / S2; + sky->elements.Flt[iy] = So; + backgnd->elements.Flt[iy] = 0.5*(bckF_lf[iy] + bckF_rt[iy]); + } + + free (tmpv); + free (tmpx); + + free (bckF_lf); + free (bckF_rt); + free (skyF_lf); + free (skyF_rt); + + free (bckX_lf); + free (bckX_rt); + free (skyX_lf); + free (skyX_rt); + + return TRUE; + + usage: + gprint (GP_ERR, "USAGE: deimos getobj (buffer) (Xref) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n"); + // can I make default values for trace (0.0), profile (?), psf (fraction of profile?), stilt (0.0) + return FALSE; +} + + /* + // convolve the object flux at this y-coord with PSF and add sky + for (int ix = 0; i < Nx; i++) { + opihi_flt g = s = 0; + for (int n = -Npsf; n <= Npsf; n++) { + if (ix + n < 0) continue; + if (ix + n >= Nx) continue; + opihi_flt value = objV ? objV[ix+n] : 1.0; + s += psfV[n]*value; + g += psfV[n]; + } + out[ix + iy*Nx] = Npsf ? s / g : (objV ? objV[ix+n] : 1.0); + if (skyV) { + out[ix + iy*Nx] += skyV[iy] + } + */ + Index: trunk/Ohana/src/opihi/mana/deimos_mkalt.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_mkalt.c (revision 41158) +++ trunk/Ohana/src/opihi/mana/deimos_mkalt.c (revision 41159) @@ -1,4 +1,6 @@ # include "data.h" # include "deimos.h" + +XXX: deprecated // Use this to test the deimos_make_object function used by deimos_fitobj @@ -92,5 +94,5 @@ // generate the model based on the test values - float *model = deimos_make_model (objV, skyV, bckV, psf, profile, trace, Nx, Ny, 0); + float *model = deimos_make_model (objV, skyV, bckV, psf, profile, trace, NULL, Nx, Ny, 0); ResetBuffer (output, Nx, Ny, -32, 0.0, 1.0); Index: trunk/Ohana/src/opihi/mana/deimos_mkmodel.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_mkmodel.c (revision 41158) +++ trunk/Ohana/src/opihi/mana/deimos_mkmodel.c (revision 41159) @@ -86,11 +86,12 @@ // generate a subimage of the full model for the range row - row + Ny // obj,sky,bck are subset vectors corresponding to the range row to row+Ny -float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *trace, int Nx, int Ny, int row) { - - float *raw = deimos_make_straight_image (obj, sky, psf, Nx, Ny, row); +// the slit_trace_* and psf_trace refer to the full image starting at row +float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, Spline *psf_trace, int Nx, int Ny, int row) { + + float *raw = deimos_make_straight_image (obj, sky, psf, psf_trace, Nx, Ny, row); float *new = deimos_apply_tilt (raw, Nx, Ny); deimos_apply_profile (profile, new, Nx, Ny); deimos_add_background (bck, new, Nx, Ny, row); - float *out = deimos_apply_trace (trace, new, Nx, Ny, row); + float *out = deimos_apply_trace (slit_trace_red, slit_trace_blu, redlimit, new, Nx, Ny, row); free (raw); @@ -104,5 +105,5 @@ // *** generate a raw image with no slit tilt and no trace offsets. apply the PSF to the // object flux, add in the sky flux -float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row) { +float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row) { OHANA_UNUSED_PARAM(row); @@ -113,5 +114,5 @@ // Xref is a global supplied by the user (defaults to Nx/2 if < 0) - int Xoff = (int)(Xref - Npsf); + int XoffRef = (int)(Xref - Npsf); // Nx : width of output window @@ -123,7 +124,26 @@ // we are generating the image for just the row range row to row+Ny + // note that obj & sky are subset vectors for just this range of pixels opihi_flt objVy = obj[iy]; // if we shift to using the row offset: obj[iy + row] opihi_flt skyVy = sky[iy]; // if we shift to using the row offset: sky[iy + row] + // integral and fractional pixel offsets of PSF (none, by default) + int dxi = 0; + float dxf = 0.0, dxr = 1.0; + int Xoff = XoffRef; + + if (psf_trace) { + // the psf_trace is referenced against the full image, so we need to add row to iy: + float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy + row); + dxi = floor(dx); + dxf = dx - dxi; // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7 + dxr = 1 - dxf; + Xoff += dxi; + // fprintf (stderr, "%d: %f : %d\n", iy + row, dx, Xoff); + } + + // if fractional offset is small, do not interpolate + int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE; + // flux = obj * PSF + sky for (int ix = 0; ix < Nx; ix++) { @@ -136,5 +156,20 @@ // only add in the flux if we are in range of the PSF if ((n >= 0) && (n < NpsfFull)) { - value += objVy * psfV[n]; + if (doInterp) { + if ((n > 0) && (n < NpsfFull)) { + if (isfinite(psfV[n]) && isfinite(psfV[n-1])) { + value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf; + } + } + if (n == 0) { + if (isfinite(psfV[n])) { + value += objVy*psfV[n]*dxr; + } + } + } else { + if (isfinite(psfV[n])) { + value += objVy * psfV[n]; + } + } } out[ix + iy*Nx] = value; @@ -219,5 +254,5 @@ // shift pixels in x based on the trace -float *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row) { +float *deimos_apply_trace (Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, float *input, int Nx, int Ny, int row) { ALLOCATE_PTR (output, float, Nx*Ny); @@ -227,5 +262,7 @@ // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center // NOTE: spline is evaluated at full dispersion coordinate - float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy + row); + int iy_full = iy + row; + Spline *slit_trace = (iy_full < redlimit) ? slit_trace_red : slit_trace_blu; + float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy_full); // extract the integer pixel offset and the fractional offset Index: trunk/Ohana/src/opihi/mana/deimos_mkobj.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_mkobj.c (revision 41158) +++ trunk/Ohana/src/opihi/mana/deimos_mkobj.c (revision 41159) @@ -6,33 +6,60 @@ // input parameters: - // * trace : spline fit of slit central x pos vs y-coord - // * profile : slit window profile (vector) - // * object : vector of object flux vs y-coord - // * sky : vector of local sky signal vs y-coord - // * background : vector of extra-slit background flux vs y-coord - // * PSF : point-spread function vector (flux normalized, x-dir) - // * LSF : sline-spread function vector (flux normalized, y-dir) - // * stilt : slit tilt response : 2D kernel? + // * slit_trace_red : spline fit of slit central x pos vs y-coord for red chip + // * slit_trace_blu : spline fit of slit central x pos vs y-coord for blu chip + // * psf_trace : spline fit of psf central x pos vs y-coord + // * profile : slit window profile (vector) + // * object : vector of object flux vs y-coord + // * sky : vector of local sky signal vs y-coord + // * backgnd : vector of extra-slit background flux vs y-coord + // * PSF : point-spread function vector (flux normalized, x-dir) + // * LSF : sline-spread function vector (flux normalized, y-dir) + // * stilt : slit tilt response : 2D kernel? + + // * redlimit : pixel break between red and blu chips + // * Nwave : number of wavelength pixels in output image int N; // if any of these are not defined, they will have assumed identity values - Spline *trace = NULL; - Vector *profile = NULL; - Vector *object = NULL; - Vector *sky = NULL; - Vector *backgnd = NULL; - - Vector *psf = NULL; - Vector *lsf = NULL; + Spline *slit_trace_red = NULL; + Spline *slit_trace_blu = NULL; + Spline *psf_trace = NULL; + Vector *profile = NULL; + Vector *object = NULL; + Vector *sky = NULL; + Vector *backgnd = NULL; + + Vector *psf = NULL; + Vector *lsf = NULL; // add tilt later - Buffer *output = NULL; + Buffer *output = NULL; float stilt = 0.0; // angle of the slit - if ((N = get_argument (argc, argv, "-trace"))) { - remove_argument (N, &argc, argv); - if ((trace = FindSpline (argv[N])) == NULL) return (FALSE); + // for a red vs blu spline, we need to specify the split point + // XXX this is REALLY ad-hoc for Deimos. not sure how to make + // this more generic (need to define the ranges somewhere) + int redlimit = 4096; + if ((N = get_argument (argc, argv, "-redlimit"))) { + remove_argument (N, &argc, argv); + redlimit = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + if ((N = get_argument (argc, argv, "-slit-trace-red"))) { + remove_argument (N, &argc, argv); + if ((slit_trace_red = FindSpline (argv[N])) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-slit-trace-blu"))) { + remove_argument (N, &argc, argv); + if ((slit_trace_blu = FindSpline (argv[N])) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-psf-trace"))) { + remove_argument (N, &argc, argv); + if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE); remove_argument (N, &argc, argv); } @@ -81,6 +108,24 @@ } + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + remove_argument (N, &argc, argv); + VERBOSE = TRUE; + } + if (argc != 3) { - gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-Nwave N] [-object vector] [-sky vector] [-backgnd vector] [-trace spline] [-profile vector]\n"); + gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-object vector] [-sky vector] [-backgnd vector] [-slit-trace-red spline] [-slit-trace-blu spline] [-psf-trace spline] [-profile vector] [-psf vector] [-lsf vector] [-stilt angle] [-Nwave Npixel] [-redlimit pixel]\n"); + gprint (GP_ERR, " slit_trace_red : spline fit of slit central x pos vs y-coord for red chip\n"); + gprint (GP_ERR, " slit_trace_blu : spline fit of slit central x pos vs y-coord for blu chip\n"); + gprint (GP_ERR, " psf_trace : spline fit of psf central x pos vs y-coord\n"); + gprint (GP_ERR, " profile : slit window profile (vector)\n"); + gprint (GP_ERR, " object : vector of object flux vs y-coord \n"); + gprint (GP_ERR, " sky : vector of local sky signal vs y-coord\n"); + gprint (GP_ERR, " backgnd : vector of extra-slit background flux vs y-coord\n"); + gprint (GP_ERR, " PSF : point-spread function vector (flux normalized, x-dir)\n"); + gprint (GP_ERR, " LSF : sline-spread function vector (flux normalized, y-dir)\n"); + gprint (GP_ERR, " stilt : slit tilt response : 2D kernel? \n"); + gprint (GP_ERR, " redlimit : pixel break between red and blu chips\n"); + gprint (GP_ERR, " Nwave : number of wavelength pixels in output image\n"); return FALSE; } @@ -123,10 +168,11 @@ // if we are appying a trace offset spline, we need to generate an output window which // is Nx + the full swing of the trace, then window back down - if (trace) { + if (slit_trace_red && slit_trace_blu) { float dXmin = +1000; float dXmax = -1000; for (int iy = 0; iy < Ny; iy++) { // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center - float dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy); + Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu; + float dx = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy); dXmin = MIN (dx, dXmin); dXmax = MAX (dx, dXmax); @@ -181,4 +227,21 @@ if (psf) Npof -= Npsf; + // integral and fractional pixel offsets of PSF (none, by default) + int dxi = 0; + float dxf = 0.0, dxr = 1.0; + + if (psf_trace) { + float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy); + dxi = floor(dx); + dxf = dx - dxi; // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7 + dxr = 1 - dxf; + Npof += dxi; + // fprintf (stderr, "%d: %f : %d\n", iy, dx, Npof); + if (VERBOSE && (iy % 200 == 0)) { gprint (GP_ERR, "%d : %f : %d : %f %f\n", iy, dx, dxi, dxf, dxr); } + } + + // if fractional offset is small, do not interpolate + int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE; + // flux = obj * PSF + sky for (int ix = 0; ix < Nx; ix++) { @@ -191,5 +254,20 @@ // only add in the flux if we are in range of the PSF if ((n >= 0) && (n < NpsfFull)) { - value += objVy * psfV[n]; + if (doInterp) { + if ((n > 0) && (n < NpsfFull)) { + if (isfinite(psfV[n]) && isfinite(psfV[n-1])) { + value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf; + } + } + if (n == 0) { + if (isfinite(psfV[n])) { + value += objVy*psfV[n]*dxr; + } + } + } else { + if (isfinite(psfV[n])) { + value += objVy*psfV[n]; + } + } } } else { @@ -329,5 +407,5 @@ // shift pixels in x based on the trace - if (trace) { + if (slit_trace_red && slit_trace_blu) { ALLOCATE_PTR (outTraceBuffer, char, NxBase*Ny*sizeof(float)); float *outTrace = (float *) outTraceBuffer; @@ -337,5 +415,6 @@ // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center - float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy); + Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu; + float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy); // extract the integer pixel offset and the fractional offset Index: trunk/Ohana/src/opihi/mana/deimos_mkslit.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_mkslit.c (revision 41158) +++ trunk/Ohana/src/opihi/mana/deimos_mkslit.c (revision 41159) @@ -7,14 +7,21 @@ // input parameters: // slit profile response : vector of fractional flux vs x-coord - // trace : spline fit of slit central x pos vs y-coord - // flux : vector of input signal vs y-coord + // trace_ref,trace_blu : spline fit of slit central x pos vs y-coord + // flux (optional) : vector of input signal vs y-coord + // sky (optional) : vector of input signal vs y-coord + // fluxbins (optional) : size of output image in y-direction + // NOTE: one of flux, sky, or fluxbins must be defined + + // output : buff (an image with Nx defined by profile and Ny defined + // by flux, sky, or fluxbins) int N; - Vector *profile = NULL; - Vector *flux = NULL; - Vector *sky = NULL; - Spline *trace = NULL; - Buffer *buff = NULL; + Vector *profile = NULL; + Vector *flux = NULL; + Vector *sky = NULL; + Spline *trace_red = NULL; + Spline *trace_blu = NULL; + Buffer *buff = NULL; // user-specified flux value (otherwise assumed to be 1) @@ -40,6 +47,24 @@ } - if (argc != 4) { - gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n"); + // either flux or fluxbins must be specified to define output size + int profilebins = 0; + if ((N = get_argument (argc, argv, "-profilebins"))) { + remove_argument (N, &argc, argv); + profilebins = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + // for a red vs blu spline, we need to specify the split point + // XXX this is REALLY ad-hoc for Deimos. not sure how to make + // this more generic (need to define the ranges somewhere) + int redlimit = 4096; + if ((N = get_argument (argc, argv, "-redlimit"))) { + remove_argument (N, &argc, argv); + redlimit = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 5) { + gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace_red) (trace_blu) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n"); return FALSE; } @@ -51,10 +76,12 @@ // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior - if ((profile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); - if ((trace = FindSpline (argv[2])) == NULL) return (FALSE); - if ((buff = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE); + if ((profile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((trace_red = FindSpline (argv[2])) == NULL) return (FALSE); + if ((trace_blu = FindSpline (argv[3])) == NULL) return (FALSE); + if ((buff = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE); // define the output window - int Nx = 2*profile[0].Nelements + 1; + // if no profile half-width is specified, use supplied profile + int Nx = profilebins ? 2*profilebins + 1 : 2*profile[0].Nelements + 1; int Ny = flux ? flux[0].Nelements : fluxbins; @@ -73,4 +100,5 @@ // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center + Spline *trace = (iy < redlimit) ? trace_red : trace_blu; float dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy); @@ -90,5 +118,5 @@ // set the sky in the entire window - out[ix + iy*Nx] = sky ? sky->elements.Flt[iy] : 0.0; + out[ix + iy*Nx] = NAN; // equivalent coord in the profile: @@ -98,4 +126,5 @@ // a default value: + float vsky = sky ? sky->elements.Flt[iy] : 0.0; float vout = NAN; @@ -121,5 +150,5 @@ if (flux) vout *= flux->elements.Flt[iy]; - out[ix + iy*Nx] += vout; + out[ix + iy*Nx] = vsky + vout; } } Index: trunk/Ohana/src/opihi/mana/test/deimos_fitarc.sh =================================================================== --- trunk/Ohana/src/opihi/mana/test/deimos_fitarc.sh (revision 41159) +++ trunk/Ohana/src/opihi/mana/test/deimos_fitarc.sh (revision 41159) @@ -0,0 +1,91 @@ + +macro test_01 + if ($0 != 2) + echo "USAGE: test_01 (offset)" + break + end + + # make a gaussian profile + create x 0 31 + set yo = 105.1*exp(-0.5*(x - 12.2)^2/2.2^2) + gaussdev dg yo[] 0 1 + # set dg = zero(yo) + set y = yo + dg*sqrt(yo) + + # fitarc version + dimenup y im 1 31 + deimos fitarc im 0 0 1 31 -resid imf -max-iterations 0 + dimendown imf yf1 + + echo $C0 $C2 $C1 $C3 + + $C0 = $C0 + $1 + vgauss -apply x yf2 + + # dy1 : model from fitarc + # dy2 : rebuilt model + # NOTE: vgauss and fitarc output results match exactly + # asymmetry is a result of the fitted values + + set dy1 = y - yf1 + set dy2 = y - yf2 + + lim -n 1 x dy2; clear; box; + + plot x dy1 -c red -x hist + plot x dy2 -c blue -x hist + + lim -n 2 x y; clear; box; + + plot x y -c black -x hist -lw 4 + + plot x yf1 -c red -x hist -lw 3 + plot x yf2 -c blue -x hist -lw 2 +end + +macro test_02 + if ($0 != 2) + echo "USAGE: test_02 (offset)" + break + end + + # make a gaussian profile + create x 0 31 + set yo = 105.1*exp(-0.5*(x - 12.2)^2/2.2^2) + gaussdev dg yo[] 0 1 + set dg = zero(yo) + set y = yo + dg*sqrt(yo) + + # vgauss version + $C0 = 15 + $C1 = 2 + $C2 = 100 + $C3 = 0 + vgauss x y con yf1 -q + echo $C0 $C2 $C1 $C3 + $C0 = $C0 + $1 + vgauss -apply x yf3 + + set dy1 = y - yf1 + set dy3 = y - yf3 + + # fitarc version + dimenup y im 1 31 + deimos fitarc im 0 0 1 31 -resid dim -max-iterations 0 + dimendown dim yf2 + set dy2 = y - yf2 + + lim -n 1 x dy2; clear; box; + + plot x dy1 -c red -x hist + plot x dy2 -c blue -x hist + plot x dy3 -c black -x hist + + lim -n 2 x y; clear; box; + + plot x y -c black -x hist -lw 3 + + plot x yf2 -c blue -x hist -lw 3 + plot x yf1 -c red -x hist -lw 2 + plot x yf3 -c blue60 -x hist +end Index: trunk/Ohana/src/opihi/mana/test/deimos_fitprofile.sh =================================================================== --- trunk/Ohana/src/opihi/mana/test/deimos_fitprofile.sh (revision 41159) +++ trunk/Ohana/src/opihi/mana/test/deimos_fitprofile.sh (revision 41159) @@ -0,0 +1,28 @@ + +macro memtest1 + # generate slit profile: + create xprof 0 17 + $Xpmin = 3 + $Xpmax = 13 + set profile_left = (xprof <= $Xpmin) ? exp(-0.5*(xprof - $Xpmin)^2/1.0^2) : zero(xprof) + set profile_right = (xprof >= $Xpmax) ? exp(-0.5*(xprof - $Xpmax)^2/1.0^2) : zero(xprof) + set profile_mid = (xprof > $Xpmin) && (xprof < $Xpmax) ? zero(xprof) + 1 : zero(xprof) + set fprofile = profile_left + profile_mid + profile_right + set wprofile = 0.02*fprofile + + $i = 0 + deimos fitprofile -q xprof fprofile wprofile $Xpmin $Xpmax + + # tv -n 0 out 250 750 + # center 15 50 20 + # resize 850 2050 + # tvcolor rainbow + # break + + memory check + for i 0 1000 + deimos fitprofile -q xprof fprofile wprofile $Xpmin $Xpmax + end + memory check +end + Index: trunk/Ohana/src/opihi/mana/test/deimos_fitslit.sh =================================================================== --- trunk/Ohana/src/opihi/mana/test/deimos_fitslit.sh (revision 41159) +++ trunk/Ohana/src/opihi/mana/test/deimos_fitslit.sh (revision 41159) @@ -0,0 +1,31 @@ + +macro memtest1 + # make a delta-function profile + create x 0 11 + set yg = exp(-0.25*(x - 5)^2) + + # flux vs wavelength is flat & short + create w 0 100 + set f = zero(w) + 1 + + delete -q Xn Yn + vlist Xn 0 100 + vlist Yn 0 2 + + spline create trace Xn Yn + spline apply trace w dX + $i = 0 + + # run once to ensure the output is create + deimos mkslit yg trace trace outbuff -flux f + mgaussdev noise outbuff[][0] outbuff[0][] 0.0 0.2 + + set obs = 2*outbuff + noise + deimos fitslit obs outbuff outf outs + + memory check + for i 0 1000 + deimos mkslit yg trace trace outbuff -flux f + end + memory check +end Index: trunk/Ohana/src/opihi/mana/test/deimos_mkalt.sh =================================================================== --- trunk/Ohana/src/opihi/mana/test/deimos_mkalt.sh (revision 41159) +++ trunk/Ohana/src/opihi/mana/test/deimos_mkalt.sh (revision 41159) @@ -0,0 +1,55 @@ + +macro memtest1 + # generate LSF: + create x 0 7 + set lsf = exp(-0.5*(x-3)^2) + + # generate PSF: + create xp 0 9 + set psf = exp(-0.5*(xp-5)^2/1.5^2) + + # generate slit profile: + create xprof 0 17 + $Xpmin = 3 + $Xpmax = 13 + set profile_left = (xprof <= $Xpmin) ? exp(-0.5*(xprof - $Xpmin)^2/1.0^2) : zero(xprof) + set profile_right = (xprof >= $Xpmax) ? exp(-0.5*(xprof - $Xpmax)^2/1.0^2) : zero(xprof) + set profile_mid = (xprof > $Xpmin) && (xprof < $Xpmax) ? zero(xprof) + 1 : zero(xprof) + set profile = profile_left + profile_mid + profile_right + + # generate object flux: + create wave 0 100 + set flux = 100 + 0.5*wave + + # generate sky flux: + set sky = 200 + zero(wave) + + # off-slit background: + set backgnd = 300 + zero(wave) + + # slit trace + vlist Xn 0 100 + vlist Yn 0 4 + spline create slit_trace_in Xn Yn + + # psf trace + vlist Xn 0 25 50 75 100 + vlist Yn 0 2 4 2 0 + spline create psf_trace_in Xn Yn + + $i = 0 + deimos mkalt out 31 -psf psf -object flux -sky sky -stilt 20.0 -profile profile -backgnd backgnd -trace slit_trace_in + + # tv -n 0 out 250 750 + # center 15 50 20 + # resize 850 2050 + # tvcolor rainbow + # break + + memory check + for i 0 1000 + deimos mkalt out 31 -psf psf -object flux -sky sky -stilt 20.0 -profile profile -backgnd backgnd -trace slit_trace_in + end + memory check +end + Index: trunk/Ohana/src/opihi/mana/test/deimos_mkobj.sh =================================================================== --- trunk/Ohana/src/opihi/mana/test/deimos_mkobj.sh (revision 41159) +++ trunk/Ohana/src/opihi/mana/test/deimos_mkobj.sh (revision 41159) @@ -0,0 +1,55 @@ + +macro memtest1 + # generate LSF: + create x 0 7 + set lsf = exp(-0.5*(x-3)^2) + + # generate PSF: + create xp 0 9 + set psf = exp(-0.5*(xp-5)^2/1.5^2) + + # generate slit profile: + create xprof 0 17 + $Xpmin = 3 + $Xpmax = 13 + set profile_left = (xprof <= $Xpmin) ? exp(-0.5*(xprof - $Xpmin)^2/1.0^2) : zero(xprof) + set profile_right = (xprof >= $Xpmax) ? exp(-0.5*(xprof - $Xpmax)^2/1.0^2) : zero(xprof) + set profile_mid = (xprof > $Xpmin) && (xprof < $Xpmax) ? zero(xprof) + 1 : zero(xprof) + set profile = profile_left + profile_mid + profile_right + + # generate object flux: + create wave 0 100 + set flux = 100 + 0.5*wave + + # generate sky flux: + set sky = 200 + zero(wave) + + # off-slit background: + set backgnd = 300 + zero(wave) + + # slit trace + vlist Xn 0 100 + vlist Yn 0 4 + spline create slit_trace_in Xn Yn + + # psf trace + vlist Xn 0 25 50 75 100 + vlist Yn 0 2 4 2 0 + spline create psf_trace_in Xn Yn + + $i = 0 + deimos mkobj out 31 -psf psf -object flux -sky sky -stilt 20.0 -lsf lsf -profile profile -backgnd backgnd -slit-trace slit_trace_in -psf-trace psf_trace_in + + # tv -n 0 out 250 750 + # center 15 50 20 + # resize 850 2050 + # tvcolor rainbow + # break + + memory check + for i 0 1000 + deimos mkobj out 31 -psf psf -object flux -sky sky -stilt 20.0 -lsf lsf -profile profile -backgnd backgnd -slit-trace slit_trace_in -psf-trace psf_trace_in + end + memory check +end + Index: trunk/Ohana/src/opihi/mana/test/deimos_mkslit.sh =================================================================== --- trunk/Ohana/src/opihi/mana/test/deimos_mkslit.sh (revision 41159) +++ trunk/Ohana/src/opihi/mana/test/deimos_mkslit.sh (revision 41159) @@ -0,0 +1,28 @@ + + +macro memtest1 + # make a delta-function profile + create x 0 11 + set yg = exp(-0.25*(x - 5)^2) + + # flux vs wavelength is flat & short + create w 0 100 + set f = zero(w) + 1 + + delete -q Xn Yn + vlist Xn 0 100 + vlist Yn 0 2 + + spline create trace Xn Yn + spline apply trace w dX + $i = 0 + + # run once to ensure the output is create + deimos mkslit yg trace trace outbuff -flux f + + memory check + for i 0 1000 + deimos mkslit yg trace trace outbuff -flux f + end + memory check +end