Index: trunk/Ohana/src/opihi/cmd.data/nnet_train.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/nnet_train.c (revision 40641) +++ trunk/Ohana/src/opihi/cmd.data/nnet_train.c (revision 40642) @@ -3,5 +3,4 @@ // these are local only float nnet_sigmoid (float value); -void nnet_sortseq (float *X, int *Y, int N); int nnet_onelayer (Nnet *nnet, int L); void nnet_reset_Nabla (Nnet *nnet); @@ -191,5 +190,5 @@ // generate a random sequence : used to select random mini batches for (int i = 0; i < Ntrial; i++) { seq[i] = i; rnd[i] = drand48(); } - nnet_sortseq (rnd, seq, Ntrial); + sort_float_index (rnd, seq, Ntrial); int Npass = Ntrial / Nmini; @@ -514,18 +513,2 @@ return 1.0 / (1.0 + exp(-value)); } - -void nnet_sortseq (float *X, int *Y, int N) { - -# define SWAPFUNC(A,B){ float ftmp; int itmp; \ - ftmp = X[A]; X[A] = X[B]; X[B] = ftmp; \ - itmp = Y[A]; Y[A] = Y[B]; Y[B] = itmp; \ - } -# define COMPARE(A,B)(X[A] < X[B]) - - OHANA_SORT (N, COMPARE, SWAPFUNC); - -# undef SWAPFUNC -# undef COMPARE - -} - Index: trunk/Ohana/src/opihi/cmd.data/sort.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/sort.c (revision 40641) +++ trunk/Ohana/src/opihi/cmd.data/sort.c (revision 40642) @@ -1,37 +1,3 @@ # include "data.h" - -void fsortindex (opihi_flt *X, int *IDX, int N) { - -# define SWAPFUNC(A,B){ opihi_flt tmp; int itmp; \ - tmp = X[A]; X[A] = X[B]; X[B] = tmp; \ - itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \ -} - -// # define COMPARE(A,B)(X[A] < X[B]) -# define COMPARE(A,B)((!isfinite(X[A]) && isfinite(X[B])) || (X[A] < X[B])) - - OHANA_SORT (N, COMPARE, SWAPFUNC); - -# undef SWAPFUNC -# undef COMPARE - -} - -void isortindex (opihi_int *X, int *IDX, int N) { - -# define SWAPFUNC(A,B){ opihi_int tmp; int itmp; \ - tmp = X[A]; X[A] = X[B]; X[B] = tmp; \ - itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \ -} - -# define COMPARE(A,B)(X[A] < X[B]) -// # define COMPARE(A,B)((!isfinite(X[A]) && isfinite(X[B])) || (X[A] < X[B])) - - OHANA_SORT (N, COMPARE, SWAPFUNC); - -# undef SWAPFUNC -# undef COMPARE - -} // XXX add an option to NOT sort, but return an index instead? @@ -78,7 +44,7 @@ /* sort key & index */ if (vec[0][0].type == OPIHI_FLT) { - fsortindex (vec[0][0].elements.Flt, index, Nval); + sort_opihi_flt_index (vec[0][0].elements.Flt, index, Nval); } else { - isortindex (vec[0][0].elements.Int, index, Nval); + sort_int_index (vec[0][0].elements.Int, index, Nval); } Index: trunk/Ohana/src/opihi/dvo/gstar.c =================================================================== --- trunk/Ohana/src/opihi/dvo/gstar.c (revision 40641) +++ trunk/Ohana/src/opihi/dvo/gstar.c (revision 40642) @@ -71,19 +71,4 @@ void printPhotcodeSequence (Average *average, SecFilt *secfilt, int entry, int type); -void sort_index_float (float *X, int *IDX, int N) { - -# define SWAPFUNC(A,B){ float tmp; int itmp; \ - tmp = X[A]; X[A] = X[B]; X[B] = tmp; \ - itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \ - } -# define COMPARE(A,B)(X[A] < X[B]) - - OHANA_SORT (N, COMPARE, SWAPFUNC); - -# undef SWAPFUNC -# undef COMPARE - -} - int gstar (int argc, char **argv) { @@ -712,5 +697,5 @@ } - sort_index_float (value, index, catalog.average[k].Nmeasure); + sort_float_index (value, index, catalog.average[k].Nmeasure); for (j = 0; j < catalog.average[k].Nmeasure; j++, m++) { Index: trunk/Ohana/src/opihi/include/data.h =================================================================== --- trunk/Ohana/src/opihi/include/data.h (revision 40641) +++ trunk/Ohana/src/opihi/include/data.h (revision 40642) @@ -249,4 +249,9 @@ int GetKapaChannelFromString (char *string); +// in sort_funcs.c: +void sort_opihi_flt_index (opihi_flt *X, int *IDX, int N); +void sort_int_index (opihi_int *X, int *IDX, int N); +void sort_float_index (float *X, int *IDX, int N); + /*** Neural Network functions (nnet_*) ***/ Index: trunk/Ohana/src/opihi/include/deimos.h =================================================================== --- trunk/Ohana/src/opihi/include/deimos.h (revision 40642) +++ trunk/Ohana/src/opihi/include/deimos.h (revision 40642) @@ -0,0 +1,20 @@ + +typedef struct { + float chisq; + opihi_flt *obj; + opihi_flt *sky; + opihi_flt *bck; +} DeimosResult; + +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); +void deimos_make_kernel (float stilt, int Nx); +void deimos_free_kernel (); +void deimos_set_cross_ref (int value, int Nx); + +// internal to make_model: +float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row); +float *deimos_apply_tilt (float *input, int Nx, int Ny); +void deimos_apply_profile (Vector *profile, float *out, int Nx, int Ny); +void deimos_add_background (opihi_flt *backgnd, float *out, int Nx, int Ny, int row); +float *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row); + Index: trunk/Ohana/src/opihi/lib.data/Makefile =================================================================== --- trunk/Ohana/src/opihi/lib.data/Makefile (revision 40641) +++ trunk/Ohana/src/opihi/lib.data/Makefile (revision 40642) @@ -29,4 +29,5 @@ $(SDIR)/nnet.$(ARCH).o \ $(SDIR)/precess.$(ARCH).o \ +$(SDIR)/sort_funcs.$(ARCH).o \ $(SDIR)/starfuncs.$(ARCH).o \ $(SDIR)/hermitian.$(ARCH).o \ Index: trunk/Ohana/src/opihi/lib.data/sort_funcs.c =================================================================== --- trunk/Ohana/src/opihi/lib.data/sort_funcs.c (revision 40642) +++ trunk/Ohana/src/opihi/lib.data/sort_funcs.c (revision 40642) @@ -0,0 +1,52 @@ +# include "data.h" + +void sort_opihi_flt_index (opihi_flt *X, int *IDX, int N) { + +# define SWAPFUNC(A,B){ opihi_flt tmp; int itmp; \ + tmp = X[A]; X[A] = X[B]; X[B] = tmp; \ + itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \ +} + +// # define COMPARE(A,B)(X[A] < X[B]) +# define COMPARE(A,B)((!isfinite(X[A]) && isfinite(X[B])) || (X[A] < X[B])) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +void sort_int_index (opihi_int *X, int *IDX, int N) { + +# define SWAPFUNC(A,B){ opihi_int tmp; int itmp; \ + tmp = X[A]; X[A] = X[B]; X[B] = tmp; \ + itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \ +} + +# define COMPARE(A,B)(X[A] < X[B]) +// # define COMPARE(A,B)((!isfinite(X[A]) && isfinite(X[B])) || (X[A] < X[B])) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +void sort_float_index (float *X, int *IDX, int N) { + +# define SWAPFUNC(A,B){ float tmp; int itmp; \ + tmp = X[A]; X[A] = X[B]; X[B] = tmp; \ + itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \ + } +// # define COMPARE(A,B)(X[A] < X[B]) +# define COMPARE(A,B)((!isfinite(X[A]) && isfinite(X[B])) || (X[A] < X[B])) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + Index: trunk/Ohana/src/opihi/mana/Makefile =================================================================== --- trunk/Ohana/src/opihi/mana/Makefile (revision 40641) +++ trunk/Ohana/src/opihi/mana/Makefile (revision 40642) @@ -24,4 +24,8 @@ $(SRC)/deimos_mkslit.$(ARCH).o \ $(SRC)/deimos_fitslit.$(ARCH).o \ +$(SRC)/deimos_getobj.$(ARCH).o \ +$(SRC)/deimos_fitobj.$(ARCH).o \ +$(SRC)/deimos_mkalt.$(ARCH).o \ +$(SRC)/deimos_mkmodel.$(ARCH).o \ $(SRC)/findrowpeaks.$(ARCH).o Index: trunk/Ohana/src/opihi/mana/deimos.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos.c (revision 40641) +++ trunk/Ohana/src/opihi/mana/deimos.c (revision 40642) @@ -4,7 +4,13 @@ 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_fitobj (int argc, char **argv); static Command deimos_commands[] = { + {1, "fitobj", deimos_fitobj, "fit for object parameters"}, + {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"}, Index: trunk/Ohana/src/opihi/mana/deimos_fitobj.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_fitobj.c (revision 40642) +++ trunk/Ohana/src/opihi/mana/deimos_fitobj.c (revision 40642) @@ -0,0 +1,368 @@ +# include "data.h" +# include "deimos.h" + +// internal functions to fitobj +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); +opihi_flt *deimos_make_test (opihi_flt *guess, float range, int row, int Nrow, int fullInput); + +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 + + 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); + } + + // 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; + } + + // Save the best 3 results for each of obj, sky, bck: + Vector *obj1, *obj2, *obj3; + if ((obj1 = SelectVector ("obj1", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (obj1, OPIHI_FLT, Ny); + if ((obj2 = SelectVector ("obj2", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (obj2, OPIHI_FLT, Ny); + if ((obj3 = SelectVector ("obj3", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (obj3, OPIHI_FLT, Ny); + Vector *sky1, *sky2, *sky3; + if ((sky1 = SelectVector ("sky1", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (sky1, OPIHI_FLT, Ny); + if ((sky2 = SelectVector ("sky2", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (sky2, OPIHI_FLT, Ny); + if ((sky3 = SelectVector ("sky3", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (sky3, OPIHI_FLT, Ny); + Vector *bck1, *bck2, *bck3; + if ((bck1 = SelectVector ("bck1", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bck1, OPIHI_FLT, Ny); + if ((bck2 = SelectVector ("bck2", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bck2, OPIHI_FLT, Ny); + if ((bck3 = SelectVector ("bck3", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bck3, OPIHI_FLT, Ny); + + for (int iy = 0; iy < Ny; iy++) { + obj1->elements.Flt[iy] = 0.0; sky1->elements.Flt[iy] = 0.0; bck1->elements.Flt[iy] = 0.0; + obj2->elements.Flt[iy] = 0.0; sky2->elements.Flt[iy] = 0.0; bck2->elements.Flt[iy] = 0.0; + obj3->elements.Flt[iy] = 0.0; sky3->elements.Flt[iy] = 0.0; bck3->elements.Flt[iy] = 0.0; + } + + // 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; + ALLOCATE_PTR (result, DeimosResult, NRESULT); + + // work on segments of Nrow at a time + // XXX worry about last segment (< Nrow rows) + for (int row = 0; row < Ny - Nrow; row += Nrow) { + + 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. + + 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); + } + + // first step, generate Niter samples (chisq + parameters) + + // 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, objRange, row, Nrow, TRUE); + opihi_flt *skytest = deimos_make_test (skyVal, skyRange, row, Nrow, TRUE); + opihi_flt *bcktest = deimos_make_test (bckVal, bckRange, 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); + + // 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, or even concept is good + objRange *= 0.75; + skyRange *= 0.75; + bckRange *= 0.75; + + // 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, objRange, row, Nrow, FALSE); + opihi_flt *skytest = deimos_make_test (result[entry].sky, skyRange, row, Nrow, FALSE); + opihi_flt *bcktest = deimos_make_test (result[entry].bck, bckRange, 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); + + // 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 + int Nv = 0; + for (int iy = 0; iy < Nrow; iy++) { + Nv = IDX[0]; obj1->elements.Flt[iy + row] = result[Nv].obj[iy]; sky1->elements.Flt[iy + row] = result[Nv].sky[iy]; bck1->elements.Flt[iy + row] = result[Nv].bck[iy]; + Nv = IDX[1]; obj2->elements.Flt[iy + row] = result[Nv].obj[iy]; sky2->elements.Flt[iy + row] = result[Nv].sky[iy]; bck2->elements.Flt[iy + row] = result[Nv].bck[iy]; + Nv = IDX[2]; obj3->elements.Flt[iy + row] = result[Nv].obj[iy]; sky3->elements.Flt[iy + row] = result[Nv].sky[iy]; bck3->elements.Flt[iy + row] = result[Nv].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 *******************/ + +opihi_flt *deimos_make_test (opihi_flt *guess, float range, 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) + + 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++) { + value[iv] = range*(drand48() - 0.5) + guess[i]; + } + return value; +} + + +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; +} Index: trunk/Ohana/src/opihi/mana/deimos_getobj.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_getobj.c (revision 40641) +++ trunk/Ohana/src/opihi/mana/deimos_getobj.c (revision 40642) @@ -1,3 +1,32 @@ # include "data.h" + +// 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; +} int deimos_getobj (int argc, char **argv) { @@ -74,7 +103,31 @@ } else { goto usage; } - if (argc != 2) 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 @@ -82,266 +135,164 @@ // based on stilt. Use these sky vectors to predict the sky under the PSF - // Here are the code stages: - - // ********* - - opihi_flt *objV = (opihi_flt *) (object ? object->elements.Flt : NULL); - opihi_flt *skyV = (opihi_flt *) (sky ? sky->elements.Flt : NULL); - opihi_flt *psfV = (opihi_flt *) (psf ? psf->elements.Flt : NULL); - opihi_flt *lsfV = (opihi_flt *) (lsf ? lsf->elements.Flt : NULL); - - // psf vector must have odd number of pixels: - int Npsf = 0; - int NpsfFull = 0; - if (psf) { - if (psf->Nelements % 2 == 0) { - gprint (GP_ERR, "psf vector must have an odd number of pixels\n"); - return FALSE; - } - NpsfFull = psf->Nelements; - Npsf = NpsfFull / 2; - } - - // lsf vector must have odd number of pixels: - int Nlsf = 0; - int NlsfFull = 0; - if (lsf) { - if (lsf->Nelements % 2 == 0) { - gprint (GP_ERR, "lsf vector must have an odd number of pixels\n"); - return FALSE; - } - NlsfFull = lsf->Nelements; - Nlsf = NlsfFull / 2; - } - - // *** apply the PSF to the object flux, add in the sky flux - - ALLOCATE_PTR (s1, float, Nx*Ny); + // *** 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 // loop over the y positions for (int iy = 0; iy < Ny; iy++) { - opihi_flt objVy = objV ? objV[iy] : 1.0; // if object is not supplied, flux defaults to 1.0 - opihi_flt skyVy = skyV ? skyV[iy] : 0.0; // if sky is not supplied, flux defaults to 0.0 + // 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 ? + } + + // 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 dXtrace = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy); - int Npof = floor(Nx/2); - if (psf) Npof -= Npsf; - - // flux = obj * PSF + sky - for (int ix = 0; ix < Nx; ix++) { - - opihi_flt value = skyVy; - - if (psfV) { - int n = ix - Npof; - // n is the pixel in the PSF corresponding to the ix pixel in the output image - // only add in the flux if we are in range of the PSF - if ((n >= 0) && (n < NpsfFull)) { - value += objVy * psfV[n]; - } - } else { - if (ix == Nx / 2) value += objVy; - } - s1[ix + iy*Nx] = value; + // 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; + + // 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; } - } - - // *** apply the slit tilt & lsf - - // first, generate the interpolation kernels. For a slit tilt of theta, there is a - // displacement in the y-direction of dy = dx sin(theta) where dx is the x-coord - // relative to the slit window center (Nx / 2). we thus need an image of size Nx, - // Nx*sin(theta) - - float sin_stilt = sin(stilt*RAD_DEG); - int Nyk = ceil(Nx * fabs(sin_stilt)); - if (Nyk < 3) Nyk = 3; // minimum of 3 pixels (or kernel assumption below fails) - Nyk += NlsfFull; - if (Nyk % 2 == 0) Nyk ++; // force Nyk to be odd - int Nyk2 = Nyk/2; // Nyk2 is the center pixel of the interpolations - ALLOCATE_PTR (kernel, float, Nx*Nyk); - int Nkpix = Nx*Nyk; - - int Nx2 = floor(Nx/2); - for (int ix = 0; ix < Nx; ix++) { - for (int iy = 0; iy < Nyk; iy++) { kernel[ix + iy*Nx] = 0.0; } - - // displacement in y-dir due to slit tilt: - float dy = (ix - Nx2) * sin_stilt; - int dyi = floor(dy); - float dyf = dy - dyi; - - if (lsf) { - int Sy = Nyk2 - Nlsf + dyi; - int doInterp = (fabs(dyf) < 1e-5) ? FALSE : TRUE; - for (int iy = 0; iy < Nyk; iy++) { - int py = iy - Sy; - if (py < 0) continue; - if (py >= NlsfFull) continue; - - float vout = NAN; - if (doInterp) { - if ((py > 0) && (py < NlsfFull - 1)) { - vout = lsfV[py]*(1 - dyf) + lsfV[py-1]*dyf; - } - if (py == 0) { - vout = lsfV[py]*dyf; - } - if (py == NlsfFull - 1) { - vout = lsfV[py]*(1.0 - dyf); - } - } else { - vout = lsfV[py]; - } - int Nkpix_i = ix + iy*Nx; - myAssert (Nkpix_i < Nkpix, "oops"); - kernel[Nkpix_i] = vout; - } - } else { - int Nkpix_i; - Nkpix_i = ix + (dyi + Nyk2 + 0)*Nx; - myAssert (Nkpix_i < Nkpix, "oops"); - kernel[Nkpix_i] = 1.0 - dyf; - - Nkpix_i = ix + (dyi + Nyk2 + 1)*Nx; - myAssert (Nkpix_i < Nkpix, "oops"); - kernel[Nkpix_i] = dyf; - fprintf (stderr, "%d, %d: %f %f\n", ix, (dyi + Nyk2), kernel[ix + (dyi + Nyk2 + 0)*Nx], kernel[ix + (dyi + Nyk2 + 1)*Nx]); + + // 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 < Npsf; 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]) } - } - - float *out = (float *) output[0].matrix.buffer; - - // next, apply the interpolation kernel to the image - - // loop over the y positions - for (int iy = 0; iy < Ny; iy++) { - - // use the tilt kernel to interpolate to this x coordinate - for (int ix = 0; ix < Nx; ix++) { - - // if the LSF is defined, we cannot bypass the convolution - if (!lsf && (fabs(stilt) < 1e-5)) { - out[ix + iy*Nx] = s1[ix + iy*Nx]; - } else { - opihi_flt g = 0.0, s = 0.0; - for (int n = -Nyk2; n <= Nyk2; n++) { - if (iy + n < 0) continue; // bottom edge of full image - if (iy + n >= Ny) continue; // top edge of full image - int Nkpix_i = ix + (n + Nyk2)*Nx; - myAssert (Nkpix_i < Nkpix, "oops"); - s += kernel[Nkpix_i]*s1[ix + (iy + n)*Nx]; - g += kernel[Nkpix_i]; - } - out[ix + iy*Nx] = s / g; - } - } - } - free (kernel); - free (s1); - - // next, apply the profile (output = input * profile) - opihi_flt *profileV = (opihi_flt *) (profile ? profile->elements.Flt : NULL); - if (profile) { - int Nprof = profile->Nelements; - int Nprof2 = Nprof / 2; - - // loop over the y positions - for (int iy = 0; iy < Ny; iy++) { - - int Sx = (int)(Nx2 - Nprof2); - - for (int ix = 0; ix < Nx; ix++) { - - // equivalent coord in the profile: - int px = ix - Sx; - if ((px >= 0) && (px < Nprof)) { - out[ix + iy*Nx] *= profileV[px]; - } else { - out[ix + iy*Nx] = 0.0; - } - } - } - } - - // add in the background - opihi_flt *backgndV = (opihi_flt *) (backgnd ? backgnd->elements.Flt : NULL); - if (backgnd) { - for (int iy = 0; iy < Ny; iy++) { - for (int ix = 0; ix < Nx; ix++) { - out[ix + iy*Nx] += backgndV[iy]; - } - } - } - - // shift pixels in x based on the trace - if (trace) { - ALLOCATE_PTR (outTraceBuffer, char, NxBase*Ny*sizeof(float)); - float *outTrace = (float *) outTraceBuffer; - int NxBase2 = NxBase / 2; - - 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); - - // extract the integer pixel offset and the fractional offset - int dxi = floor(dx); // -1.7 -> -2, -0.5 -> -1, +0.5 -> 0, +1.7 -> 1 - float dxf = dx - dxi; // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7 - float dxr = 1 - dxf; - - // if fractional offset is small, do not interpolate - int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE; - - // How do I handle this? define a temporary window which is a fraction of the full output window? - int Sx = Nx2 - NxBase2 + dxi; - - // save the original output row - // for (int ix = 0; ix < Nx; ix++) { tmprow[ix] = out[ix + iy*Nx]; } - - for (int ix = 0; ix < NxBase; ix++) { - - // equivalent coord in temporary buffer (Nx * Ny): - int px = ix + Sx; - if (px < 0) { - outTrace[ix + iy*NxBase] = out[iy*Nx]; - continue; - } - if (px >= Nx) { - outTrace[ix + iy*NxBase] = out[(Nx - 1) + iy*Nx]; - continue; - } - - // a default value: - float vout = NAN; - - int Npix = px + iy*Nx; - - if (doInterp) { - if ((px > 0) && (px < Nx - 1)) { - vout = out[Npix]*dxr + out[Npix + 1]*dxf; - } - if (px == 0) { - vout = out[Npix]; - } - if (px == Nx - 1) { - vout = out[Npix]; - } - } else { - vout = out[Npix]; - } - outTrace[ix + iy*NxBase] = vout; - } - } - - ResetBuffer (output, NxBase, Ny, -32, 0.0, 1.0); - free (output[0].matrix.buffer); - output[0].matrix.buffer = outTraceBuffer; - } + + 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; usage: - gprint (GP_ERR, "USAGE: deimos getobj (buffer) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n"); + 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; Index: trunk/Ohana/src/opihi/mana/deimos_mkalt.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_mkalt.c (revision 40642) +++ trunk/Ohana/src/opihi/mana/deimos_mkalt.c (revision 40642) @@ -0,0 +1,107 @@ +# include "data.h" +# include "deimos.h" + +// Use this to test the deimos_make_object function used by deimos_fitobj + +int deimos_mkalt (int argc, char **argv) { + + // generate model observed flux for an object in a slit with background + + // 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) + // * stilt : slit tilt response : 2D kernel? + + int N; + + // if any of these are not defined, they will have assumed identity values + Vector *obj = NULL; + Vector *sky = NULL; + Vector *bck = NULL; + + Spline *trace = NULL; + Vector *profile = NULL; + Vector *psf = NULL; + float stilt = 0.0; // angle of the slit + + Buffer *output = NULL; + + 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, "-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 ((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; } + + if (argc != 3) goto usage; + + int Nx = atoi (argv[2]); + + // all supplied vectors must be consistent + int Ny = obj->Nelements; + 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; + } + + if ((output = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE); + + opihi_flt *objV = obj->elements.Flt; + opihi_flt *skyV = sky->elements.Flt; + opihi_flt *bckV = bck->elements.Flt; + + deimos_make_kernel (stilt, Nx); + deimos_set_cross_ref (-1, Nx); + + // generate the model based on the test values + float *model = deimos_make_model (objV, skyV, bckV, psf, profile, trace, Nx, Ny, 0); + + ResetBuffer (output, Nx, Ny, -32, 0.0, 1.0); + free (output[0].matrix.buffer); + output[0].matrix.buffer = (char *) model; + + deimos_free_kernel(); + + return TRUE; + +usage: + gprint (GP_ERR, "USAGE: deimos mkalt (buffer) (Nx) [-object vector] [-sky vector] [-backgnd vector] [-trace spline] [-profile vector] [-psf psf] [-stilt tilt]\n"); + return FALSE; +} Index: trunk/Ohana/src/opihi/mana/deimos_mkmodel.c =================================================================== --- trunk/Ohana/src/opihi/mana/deimos_mkmodel.c (revision 40642) +++ trunk/Ohana/src/opihi/mana/deimos_mkmodel.c (revision 40642) @@ -0,0 +1,280 @@ +# include "data.h" +# include "deimos.h" + +// global parameters of the slit kernel: +static int Nkernel_pix = 0; // total number of pixels +static int Nkernel_y2 = 0; // center of the interpolation +static float *kernel = NULL; +static int SLIT_NO_INTERP = FALSE; // if slit tilt is minimal, do not interpolate + +static int Xref = -1; // center of the slit profile + +/****************** make_object & support functions *******************/ + +void deimos_set_cross_ref (int value, int Nx) { + Xref = value; + if (Xref < 0) { Xref = floor(Nx/2); } +} + +// kernel for slit tilt; no LSF for now +void deimos_make_kernel (float stilt, int Nx) { + + // first, generate the interpolation kernels. For a slit tilt of theta, there is a + // displacement in the y-direction of dy = dx sin(theta) where dx is the x-coord + // relative to the slit window center (Nx / 2). we thus need an image of size Nx, + // Nx*sin(theta) + + SLIT_NO_INTERP = (fabs(stilt) < 1e-5); + + float sin_stilt = sin(stilt*RAD_DEG); + int Ny = ceil(Nx * fabs(sin_stilt)) + 1; // one extra pixel buffer + if (Ny < 3) Ny = 3; // minimum of 3 pixels (or kernel assumption below fails) + + if (Ny % 2 == 0) Ny ++; // force Nyk to be odd + + Nkernel_y2 = Ny/2; // Nkernel_y2 is the center pixel in the interpolation direction + + ALLOCATE (kernel, float, Nx*Ny); + Nkernel_pix = Nx*Ny; + fprintf (stderr, "allocating %d pixels (%d,%d)\n", Nkernel_pix, Nx, Ny); + + // XXX use Xref not Nx/2? + int Nx2 = floor(Nx/2); + + for (int ix = 0; ix < Nx; ix++) { + for (int iy = 0; iy < Ny; iy++) { kernel[ix + iy*Nx] = 0.0; } + + // displacement in y-dir due to slit tilt: + float dy = (ix - Nx2) * sin_stilt; + int dyi = floor(dy); + float dyf = dy - dyi; + + // without LSF, slit tilt only affects two neighbor pixels + int pix; + pix = ix + (dyi + Nkernel_y2 + 0)*Nx; + fprintf (stderr, "%d, %d, %d -- ", ix, (dyi + Nkernel_y2 + 0), pix); + myAssert (pix < Nkernel_pix, "oops (make_kernel 1)"); + kernel[pix] = 1.0 - dyf; + + pix = ix + (dyi + Nkernel_y2 + 1)*Nx; + fprintf (stderr, "%d, %d, %d -- ", ix, (dyi + Nkernel_y2 + 1), pix); + myAssert (pix < Nkernel_pix, "oops (make_kernel 2)"); + kernel[pix] = dyf; + fprintf (stderr, ": %f %f\n", kernel[ix + (dyi + Nkernel_y2 + 0)*Nx], kernel[ix + (dyi + Nkernel_y2 + 1)*Nx]); + } +} + +void deimos_free_kernel () { + free (kernel); +} + +/* deimos_make_model generates a model image for a subset of the + wavelength range, starting at 'row' and going for Ny rows. + + we have a choice here. we could: + + 1) assume obj,sky,bck are only segments corresponding to the same wavelength range + (we still need to use 'row' to get the right trace values) + + 2) require obj,sky,bck to correspond to the full wavelength range (pass in 'row' to get + the correct values). + + I think (1) simplified the situation the most. + + */ + +// generate a subimage of the full model for the range row - row + Ny +// obj,sky,bck are vectors of the full length of the model (ie, not sub-vectors) +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); + 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); + + free (raw); + free (new); + + return out; +} + +// *************** make_model support functions ***************** + +// *** 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) { + OHANA_UNUSED_PARAM(row); + + opihi_flt *psfV = psf->elements.Flt; + + int NpsfFull = psf->Nelements; + int Npsf = NpsfFull / 2; + + // Xref is a global supplied by the user (defaults to Nx/2 if < 0) + int Xoff = (int)(Xref - Npsf); + + // Nx : width of output window + // Ny : number of output rows + ALLOCATE_PTR (out, float, Nx*Ny); + + // loop over the y positions in the output image: + for (int iy = 0; iy < Ny; iy++) { + + // we are generating the image for just the row range row to row+Ny + 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] + + // flux = obj * PSF + sky + for (int ix = 0; ix < Nx; ix++) { + + opihi_flt value = skyVy; + + int n = ix - Xoff; + + // n is the pixel in the PSF corresponding to the ix pixel in the output image + // only add in the flux if we are in range of the PSF + if ((n >= 0) && (n < NpsfFull)) { + value += objVy * psfV[n]; + } + out[ix + iy*Nx] = value; + } + } + return out; +} + +// the slit tilt shifts the flux in the dispersion direction an amount which depends on +// the cross-dispertion position. the shift is not a function of row position. +float *deimos_apply_tilt (float *input, int Nx, int Ny) { + + // Nx : width of output window + // Ny : number of output rows + ALLOCATE_PTR (output, float, Nx*Ny); + + // loop over the y positions + for (int iy = 0; iy < Ny; iy++) { + + // use the tilt kernel to interpolate to this x coordinate + for (int ix = 0; ix < Nx; ix++) { + + // if the LSF is defined, we cannot bypass the convolution + if (SLIT_NO_INTERP) { + output[ix + iy*Nx] = input[ix + iy*Nx]; + } else { + opihi_flt g = 0.0, s = 0.0; + for (int n = -Nkernel_y2; n <= Nkernel_y2; n++) { + if (iy + n < 0) continue; // bottom edge of full image + if (iy + n >= Ny) continue; // top edge of full image + int pix = ix + (n + Nkernel_y2)*Nx; + myAssert (pix < Nkernel_pix, "oops (apply_tilt)"); + s += kernel[pix]*input[ix + (iy + n)*Nx]; + g += kernel[pix]; // optimization: save g[ix], all should be the same + } + output[ix + iy*Nx] = s / g; + } + } + } + return output; +} + +// the slit window profile modifies the flux (essentially a vignetting term). +// the profile is constant vs row (do not need a row value) +void deimos_apply_profile (Vector *profile, float *out, int Nx, int Ny) { + + // next, apply the profile (output = input * profile) + opihi_flt *profileV = profile->elements.Flt; + + int Nprof = profile->Nelements; + int Nprof2 = Nprof / 2; + + // Xref is a global supplied by the user (defaults to Nx/2 if < 0) + int Xoff = (int)(Xref - Nprof2); + + // loop over the y positions + for (int iy = 0; iy < Ny; iy++) { + + for (int ix = 0; ix < Nx; ix++) { + + // equivalent coord in the profile: + int pix = ix - Xoff; + if ((pix >= 0) && (pix < Nprof)) { + out[ix + iy*Nx] *= profileV[pix]; + } else { + out[ix + iy*Nx] = 0.0; + } + } + } +} + +// supplied background model covers full dispersion range, not the sub-image +void deimos_add_background (opihi_flt *bck, float *out, int Nx, int Ny, int row) { + OHANA_UNUSED_PARAM(row); + + for (int iy = 0; iy < Ny; iy++) { + for (int ix = 0; ix < Nx; ix++) { + out[ix + iy*Nx] += bck[iy]; // if we shift to using the row offset: bck[iy + row] + } + } +} + +// shift pixels in x based on the trace +float *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row) { + + ALLOCATE_PTR (output, float, Nx*Ny); + + 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 + // NOTE: spline is evaluated at full dispersion coordinate + float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy + row); + + // extract the integer pixel offset and the fractional offset + int dxi = floor(dx); // -1.7 -> -2, -0.5 -> -1, +0.5 -> 0, +1.7 -> 1 + float dxf = dx - dxi; // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7 + float dxr = 1 - dxf; + + // if fractional offset is small, do not interpolate + int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE; + + // How do I handle this? define a temporary window which is a fraction of the full output window? + // if the output buffer is not the same width as the input buffer, use this adjustment: + // Nx2 = int(Nx/2), NxBase2 = int(NxBase/2), where Nx is input width, NxBase is output width + // int Sx = Nx2 - NxBase2 + dxi; + + for (int ix = 0; ix < Nx; ix++) { + + // equivalent coord in temporary buffer (Nx * Ny): + int pix = ix + dxi; + if (pix < 0) { + output[ix + iy*Nx] = input[iy*Nx]; + continue; + } + if (pix >= Nx) { + output[ix + iy*Nx] = input[(Nx - 1) + iy*Nx]; + continue; + } + + // a default value: + float vout = NAN; + + int Npix = pix + iy*Nx; + + if (doInterp) { + if ((pix > 0) && (pix < Nx - 1)) { + vout = input[Npix]*dxr + input[Npix + 1]*dxf; + } + if (pix == 0) { + vout = input[Npix]; + } + if (pix == Nx - 1) { + vout = input[Npix]; + } + } else { + vout = input[Npix]; + } + output[ix + iy*Nx] = vout; + } + } + return output; +} +