IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41159 for trunk


Ignore:
Timestamp:
Nov 27, 2019, 11:04:51 AM (7 years ago)
Author:
eugene
Message:

improvments to deimos extraction tools

Location:
trunk/Ohana/src/opihi/mana
Files:
1 deleted
8 edited
7 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/mana/Makefile

    r40655 r41159  
    2626$(SRC)/deimos_getobj.$(ARCH).o \
    2727$(SRC)/deimos_fitobj.$(ARCH).o \
    28 $(SRC)/deimos_fitalt.$(ARCH).o \
    29 $(SRC)/deimos_mkalt.$(ARCH).o \
    3028$(SRC)/deimos_mkmodel.$(ARCH).o \
    3129$(SRC)/deimos_arclines.$(ARCH).o \
     30$(SRC)/deimos_fitarc.$(ARCH).o \
     31$(SRC)/deimos_fitprofile.$(ARCH).o \
    3232$(SRC)/findrowpeaks.$(ARCH).o
    3333

  • trunk/Ohana/src/opihi/mana/deimos.c

    r40655 r41159  
    44int deimos_fitslit (int argc, char **argv);
    55int deimos_mkobj (int argc, char **argv);
    6 int deimos_mkalt (int argc, char **argv);
    76int deimos_getobj (int argc, char **argv);
     7// int deimos_getalt (int argc, char **argv);
    88int deimos_fitobj (int argc, char **argv);
    9 int deimos_fitalt (int argc, char **argv);
     9// int deimos_fitalt (int argc, char **argv);
    1010int deimos_arclines (int argc, char **argv);
     11int deimos_fitarc (int argc, char **argv);
     12int deimos_fitprofile (int argc, char **argv);
    1113
    1214static Command deimos_commands[] = {
    13   {1, "fitobj", deimos_fitobj, "fit for object parameters"},
    14   {1, "fitalt", deimos_fitalt, "fit for object parameters using LMM"},
    15   {1, "getobj", deimos_getobj, "determine crude object parameters"},
    16   {1, "mkobj", deimos_mkobj, "make a full object image"},
    17   {1, "mkalt", deimos_mkalt, "make a full object image (uses deimos_make_object, for a test)"},
    18   {1, "mkslit", deimos_mkslit, "make a slit image"},
    19   {1, "fitslit", deimos_fitslit, "fit slit image to observed slit flux"},
    20   {1, "arclines", deimos_arclines, "detect arclines using LSF and STILT"},
     15  {1, "fitobj",     deimos_fitobj,     "fit for object parameters using LMM"},
     16  {1, "fitalt",     deimos_fitobj,     "fit for object parameters using LMM"},
     17  {1, "getobj",     deimos_getobj,     "determine crude object parameters"},
     18  {1, "getalt",     deimos_getobj,     "determine crude object parameters"},
     19  {1, "mkobj",      deimos_mkobj,      "make a full object image"},
     20  {1, "mkslit",     deimos_mkslit,     "make a slit image"},
     21  {1, "fitslit",    deimos_fitslit,    "fit slit image to observed slit flux"},
     22  {1, "arclines",   deimos_arclines,   "detect arclines using LSF and STILT"},
     23  {1, "fitarc",     deimos_fitarc,     "fit arclamp lines"},
     24  {1, "fitprofile", deimos_fitprofile, "fit slit profile"},
    2125};
    2226

  • trunk/Ohana/src/opihi/mana/deimos_arclines.c

    r40655 r41159  
    5353  deimos_make_LSF_kernel (LSF, stilt, Nx);
    5454
     55  // FILE *Fout = NULL;
     56
    5557  // loop over the rows
    5658  for (int iy = 0; iy < Ny; iy++) {
    5759    float Fsum = 0;
     60    float Ksum = 0;
     61
     62    // if (iy == 5802) Fout = fopen ("test.5802.dat", "w");
     63    // if (iy == 5815) Fout = fopen ("test.5815.dat", "w");
     64
    5865    // cross-correlation of buffer values and kernel values
    5966    for (int ix = 0; ix < Nx; ix++) {
    … …  
    6269        if ((iy + ky) < 0) continue;
    6370        if ((iy + ky) >= Ny) continue;
     71        if (!isfinite(Fin[ix + (iy + ky)*Nx])) continue;
     72        if (!isfinite(kernel[ix + ko*Nx])) continue;
    6473        Fsum += Fin[ix + (iy + ky)*Nx] * kernel[ix + ko*Nx];
     74        Ksum += kernel[ix + ko*Nx];
     75        // if ((iy == 5815) || (iy == 5802)) {
     76        //   fprintf (Fout, "%d %d %d : %f %f\n", iy, ix, ky, Fin[ix + (iy + ky)*Nx], kernel[ix + ko*Nx]);
     77        // }
    6578      }
    6679    }
     80
     81    // if ((iy == 5815) || (iy == 5802)) fclose (Fout);
     82
    6783    coord->elements.Flt[iy] = iy;
    6884    flux->elements.Flt[iy] = Fsum;

  • trunk/Ohana/src/opihi/mana/deimos_fitslit.c

    r40659 r41159  
    11# include "data.h"
     2# define IRLS_TOLERANCE 1e-4
     3
     4float weight_cauchy_square_flt (float x2);
     5static void fitflux (float *Fwind, float *Fprof, float *weight, int Nx, double *flux, double *sky);
    26
    37int deimos_fitslit (int argc, char **argv) {
    48
    5   // fitslit (window) (slit) (flux) (sky)
    6 
    7   // int N;
     9  // fitslit (window) (vari) (slit) (flux) (sky)
    810
    911  Vector *flux    = NULL;
    … …  
    1113
    1214  Buffer *wind    = NULL;
     15  Buffer *vari    = NULL;
    1316  Buffer *slit    = NULL;
    1417
    15   if (argc != 5) {
    16     gprint (GP_ERR, "USAGE: deimos fitslit (window) (slit) (flux) (sky)\n");
    17     gprint (GP_ERR, "  inputs:  window (observed 2D flux), slit (model 2D flux)\n");
     18  int N;
     19  int Niter = 10;
     20  if ((N = get_argument (argc, argv, "-irls-iter"))) {
     21    remove_argument (N, &argc, argv);
     22    Niter = atoi (argv[N]);
     23    remove_argument (N, &argc, argv);
     24  }
     25
     26  Buffer *model = NULL;
     27  if ((N = get_argument (argc, argv, "-model"))) {
     28    remove_argument (N, &argc, argv);
     29    if ((model = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE);
     30    remove_argument (N, &argc, argv);
     31  }
     32
     33  if (argc != 6) {
     34    gprint (GP_ERR, "USAGE: deimos fitslit (window) (variance) (slit) (flux) (sky)\n");
     35    gprint (GP_ERR, "  inputs:  window (observed 2D flux), variance (on 2D flux), slit (model 2D flux)\n");
    1836    gprint (GP_ERR, "  outputs: flux (best-fit 1D flux), sky (best-fit 1D background)\n");
    1937    return FALSE;
    … …  
    2240  // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
    2341  if ((wind = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    24   if ((slit = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    25   if ((flux = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
    26   if ((sky  = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
    27 
    28   // XXX confirm slit and wind have same dimensions
     42  if ((vari = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     43  if ((slit = SelectBuffer (argv[3], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     44  if ((flux = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
     45  if ((sky  = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
    2946
    3047  // define the output window
    … …  
    3249  int Ny = wind[0].matrix.Naxis[1];
    3350 
     51  // confirm wind, vari, slit have same dimensions
    3452  if (Nx != slit[0].matrix.Naxis[0]) { gprint (GP_ERR, "size mismatch in fitslit for window and slit\n"); return FALSE; }
    3553  if (Ny != slit[0].matrix.Naxis[1]) { gprint (GP_ERR, "size mismatch in fitslit for window and slit\n"); return FALSE; }
     54  if (Nx != vari[0].matrix.Naxis[0]) { gprint (GP_ERR, "size mismatch in fitslit for window and variance\n"); return FALSE; }
     55  if (Ny != vari[0].matrix.Naxis[1]) { gprint (GP_ERR, "size mismatch in fitslit for window and variance\n"); return FALSE; }
    3656
    3757  ResetVector (flux, OPIHI_FLT, Ny);
    3858  ResetVector (sky,  OPIHI_FLT, Ny);
    3959
     60  if (model) {
     61    if (!CreateBuffer (model, Nx, Ny, -32, 1.0, 0.0)) { gprint (GP_ERR, "error generating output model buffer\n"); return FALSE; }
     62  }
     63
    4064  float *Fwind = (float *) wind[0].matrix.buffer;
     65  float *Fvari = (float *) vari[0].matrix.buffer;
    4166  float *Fprof = (float *) slit[0].matrix.buffer;
     67  float *FmodOut = model ? (float *) model[0].matrix.buffer : NULL;
     68
     69  ALLOCATE_PTR (weight, float, Nx);
     70  ALLOCATE_PTR (rawwgt, float, Nx);
    4271
    4372  // loop over the rows
    4473  for (int iy = 0; iy < Ny; iy++) {
    4574
    46     // calculate elements of the chi-square for this row:
     75    double F = 0.0, S = 0.0;
    4776
    48     double R = 0.0, P = 0.0, P2 = 0.0, f = 0.0, fp = 0.0;
    49     double wt = 1.0;
    50 
     77    // set weight based on the variance
    5178    for (int ix = 0; ix < Nx; ix++) {
    52      
    53       float fxy = Fwind[ix + iy*Nx];
    54       float pxy = Fprof[ix + iy*Nx];
    55 
    56       R  += wt;
    57       P  += pxy*wt;
    58       P2 += pxy*pxy*wt;
    59       f  += fxy*wt;
    60       fp += fxy*pxy*wt;
     79      float vari = Fvari[ix + iy*Nx];
     80      weight[ix] = !isfinite(vari) || (fabs(vari) < 1e-6) ? 1.0 : 1.0 / vari;
     81      rawwgt[ix] = weight[ix];
    6182    }
    6283
    63     double det = R*P2 - P*P;
    64     double S = (P2*f - P*fp) / det;
    65     double F = (R*fp - P*f) / det;
     84    // calculate elements of the chi-square for this row:
     85    fitflux (&Fwind[iy*Nx], &Fprof[iy*Nx], weight, Nx, &F, &S);
     86
     87    // do an IRLS loop here to reject outliers
     88    int converged = FALSE;
     89    for (int iter = 0; (iter < Niter) && !converged; iter++) {
     90
     91      double Flast = F;
     92      double Slast = S;
     93
     94      // calculate weight modification based on distances (squared).
     95      // use modifier to calculate new weighted mean
     96      for (int ix = 0; ix < Nx; ix++) {
     97
     98        // F & S are model parameters; compare the model and observed values for each pixel
     99        float Fmodel = S + F*Fprof[ix + iy*Nx];
     100        float dV = (Fwind[ix + iy*Nx] - Fmodel);
     101        float d2 = SQ(dV) * rawwgt[ix];
     102       
     103        float Mod = weight_cauchy_square_flt (d2);
     104        weight[ix] = Mod * rawwgt[ix];
     105      }
     106      fitflux (&Fwind[iy*Nx], &Fprof[iy*Nx], weight, Nx, &F, &S);
     107
     108      float dF = fabs(F - Flast);
     109      float dS = fabs(S - Slast);
     110
     111      if ((dF < F * IRLS_TOLERANCE) && (dS < S * IRLS_TOLERANCE)) converged = TRUE;
     112    }
     113
     114    if (FmodOut) {
     115      for (int ix = 0; ix < Nx; ix++) {
     116        // F & S are model parameters; compare the model and observed values for each pixel
     117        float Fmodel = S + F*Fprof[ix + iy*Nx];
     118        float dV = (Fwind[ix + iy*Nx] - Fmodel);
     119        float d2 = SQ(dV) * rawwgt[ix];
     120        float Mod = weight_cauchy_square_flt (d2);
     121        if (Mod < 0.1) {
     122          FmodOut[ix + iy*Nx] = NAN;
     123        } else {
     124          FmodOut[ix + iy*Nx] = Fmodel;
     125        }
     126      }
     127    }
    66128
    67129    flux->elements.Flt[iy] = F;
    68130    sky->elements.Flt[iy] = S;
    69131  }
    70 
    71132  return TRUE;
    72133}
    73134
     135// weight is 1 / variance
     136// Fwind, Fprof are pointers to the start of a single row
     137// Fwind is the observed flux in the window
     138// Fprof is the slit profile
     139void fitflux (float *Fwind, float *Fprof, float *weight, int Nx, double *flux, double *sky) {
     140
     141  // calculate elements of the chi-square fit for this row:
     142
     143  double R = 0.0, P = 0.0, P2 = 0.0, f = 0.0, fp = 0.0;
     144
     145  for (int ix = 0; ix < Nx; ix++) {
     146     
     147    double fxy = Fwind[ix];
     148    double pxy = Fprof[ix];
     149
     150    // skip NAN / Inf pixels in either window or slit
     151    if (!isfinite(fxy)) continue;
     152    if (!isfinite(pxy)) continue;
     153
     154    double wt = (weight == NULL) ? 1.0 : weight[ix];
     155    R  += wt;
     156    P  += pxy*wt;
     157    P2 += pxy*pxy*wt;
     158    f  += fxy*wt;
     159    fp += fxy*pxy*wt;
     160  }
     161
     162  double det = R*P2 - P*P;
     163  *sky = (P2*f - P*fp) / det;
     164  *flux = (R*fp - P*f) / det;
     165}

  • trunk/Ohana/src/opihi/mana/deimos_mkalt.c

    r40642 r41159  
    11# include "data.h"
    22# include "deimos.h"
     3
     4XXX: deprecated
    35
    46// Use this to test the deimos_make_object function used by deimos_fitobj
    … …  
    9294
    9395  // generate the model based on the test values
    94   float *model = deimos_make_model (objV, skyV, bckV, psf, profile, trace, Nx, Ny, 0);
     96  float *model = deimos_make_model (objV, skyV, bckV, psf, profile, trace, NULL, Nx, Ny, 0);
    9597
    9698  ResetBuffer (output, Nx, Ny, -32, 0.0, 1.0);

  • trunk/Ohana/src/opihi/mana/deimos_mkmodel.c

    r40659 r41159  
    8686// generate a subimage of the full model for the range row - row + Ny
    8787// obj,sky,bck are subset vectors corresponding to the range row to row+Ny
    88 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) {
    89 
    90   float *raw = deimos_make_straight_image (obj, sky, psf, Nx, Ny, row);
     88// the slit_trace_* and psf_trace refer to the full image starting at row
     89float *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) {
     90
     91  float *raw = deimos_make_straight_image (obj, sky, psf, psf_trace, Nx, Ny, row);
    9192  float *new = deimos_apply_tilt (raw, Nx, Ny);
    9293  deimos_apply_profile (profile, new, Nx, Ny);
    9394  deimos_add_background (bck, new, Nx, Ny, row);
    94   float *out = deimos_apply_trace (trace, new, Nx, Ny, row);
     95  float *out = deimos_apply_trace (slit_trace_red, slit_trace_blu, redlimit, new, Nx, Ny, row);
    9596
    9697  free (raw);
    … …  
    104105// *** generate a raw image with no slit tilt and no trace offsets.  apply the PSF to the
    105106// object flux, add in the sky flux
    106 float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row) {
     107float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row) {
    107108  OHANA_UNUSED_PARAM(row);
    108109 
    … …  
    113114
    114115  // Xref is a global supplied by the user (defaults to Nx/2 if < 0)
    115   int Xoff = (int)(Xref - Npsf);
     116  int XoffRef = (int)(Xref - Npsf);
    116117
    117118  // Nx : width of output window
    … …  
    123124   
    124125    // we are generating the image for just the row range row to row+Ny
     126    // note that obj & sky are subset vectors for just this range of pixels
    125127    opihi_flt objVy = obj[iy]; // if we shift to using the row offset: obj[iy + row]
    126128    opihi_flt skyVy = sky[iy]; // if we shift to using the row offset: sky[iy + row]
    127129     
     130    // integral and fractional pixel offsets of PSF (none, by default)
     131    int dxi = 0;
     132    float dxf = 0.0, dxr = 1.0;
     133    int Xoff = XoffRef;
     134   
     135    if (psf_trace) {
     136      // the psf_trace is referenced against the full image, so we need to add row to iy:
     137      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy + row);
     138      dxi = floor(dx);
     139      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
     140      dxr = 1 - dxf;
     141      Xoff += dxi;
     142      // fprintf (stderr, "%d: %f : %d\n", iy + row, dx, Xoff);
     143    }
     144
     145    // if fractional offset is small, do not interpolate
     146    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
     147
    128148    // flux = obj * PSF + sky
    129149    for (int ix = 0; ix < Nx; ix++) {
    … …  
    136156      // only add in the flux if we are in range of the PSF
    137157      if ((n >= 0) && (n < NpsfFull)) {
    138         value += objVy * psfV[n];
     158        if (doInterp) {
     159          if ((n > 0) && (n < NpsfFull)) {
     160            if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
     161              value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
     162            }
     163          }
     164          if (n == 0) {
     165            if (isfinite(psfV[n])) {
     166              value += objVy*psfV[n]*dxr;
     167            }
     168          }
     169        } else {
     170          if (isfinite(psfV[n])) {
     171            value += objVy * psfV[n];
     172          }
     173        }
    139174      }
    140175      out[ix + iy*Nx] = value;
    … …  
    219254
    220255// shift pixels in x based on the trace
    221 float *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row) {
     256float *deimos_apply_trace (Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, float *input, int Nx, int Ny, int row) {
    222257
    223258  ALLOCATE_PTR (output, float, Nx*Ny);
    … …  
    227262    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    228263    // NOTE: spline is evaluated at full dispersion coordinate
    229     float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy + row);
     264    int iy_full = iy + row;
     265    Spline *slit_trace = (iy_full < redlimit) ? slit_trace_red : slit_trace_blu;
     266    float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy_full);
    230267     
    231268    // extract the integer pixel offset and the fractional offset

  • trunk/Ohana/src/opihi/mana/deimos_mkobj.c

    r40659 r41159  
    66
    77  // input parameters:
    8   // * trace       : spline fit of slit central x pos vs y-coord
    9   // * profile     : slit window profile (vector)
    10   // * object      : vector of object flux vs y-coord
    11   // * sky         : vector of local sky signal vs y-coord
    12   // * background  : vector of extra-slit background flux vs y-coord
    13   // * PSF         : point-spread function vector (flux normalized, x-dir)
    14   // * LSF         : sline-spread function vector (flux normalized, y-dir)
    15   // * stilt       : slit tilt response : 2D kernel?
     8  // * slit_trace_red : spline fit of slit central x pos vs y-coord for red chip
     9  // * slit_trace_blu : spline fit of slit central x pos vs y-coord for blu chip
     10  // * psf_trace      : spline fit of psf central x pos vs y-coord
     11  // * profile        : slit window profile (vector)
     12  // * object         : vector of object flux vs y-coord
     13  // * sky            : vector of local sky signal vs y-coord
     14  // * backgnd        : vector of extra-slit background flux vs y-coord
     15  // * PSF            : point-spread function vector (flux normalized, x-dir)
     16  // * LSF            : sline-spread function vector (flux normalized, y-dir)
     17  // * stilt          : slit tilt response : 2D kernel?
     18
     19  // * redlimit       : pixel break between red and blu chips
     20  // * Nwave          : number of wavelength pixels in output image
    1621
    1722  int N;
    1823
    1924  // if any of these are not defined, they will have assumed identity values
    20   Spline *trace   = NULL;
    21   Vector *profile = NULL;
    22   Vector *object  = NULL;
    23   Vector *sky     = NULL;
    24   Vector *backgnd = NULL;
    25 
    26   Vector *psf     = NULL;
    27   Vector *lsf     = NULL;
     25  Spline *slit_trace_red = NULL;
     26  Spline *slit_trace_blu = NULL;
     27  Spline *psf_trace  = NULL;
     28  Vector *profile    = NULL;
     29  Vector *object     = NULL;
     30  Vector *sky        = NULL;
     31  Vector *backgnd    = NULL;
     32                     
     33  Vector *psf        = NULL;
     34  Vector *lsf        = NULL;
    2835  // add tilt later
    2936
    30   Buffer *output  = NULL;
     37  Buffer *output     = NULL;
    3138
    3239  float stilt = 0.0; // angle of the slit
    3340
    34   if ((N = get_argument (argc, argv, "-trace"))) {
    35     remove_argument (N, &argc, argv);
    36     if ((trace = FindSpline (argv[N])) == NULL) return (FALSE);
     41  // for a red vs blu spline, we need to specify the split point
     42  // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
     43  // this more generic (need to define the ranges somewhere)
     44  int redlimit = 4096;
     45  if ((N = get_argument (argc, argv, "-redlimit"))) {
     46    remove_argument (N, &argc, argv);
     47    redlimit = atoi (argv[N]);
     48    remove_argument (N, &argc, argv);
     49  }
     50
     51  if ((N = get_argument (argc, argv, "-slit-trace-red"))) {
     52    remove_argument (N, &argc, argv);
     53    if ((slit_trace_red = FindSpline (argv[N])) == NULL) return (FALSE);
     54    remove_argument (N, &argc, argv);
     55  }
     56  if ((N = get_argument (argc, argv, "-slit-trace-blu"))) {
     57    remove_argument (N, &argc, argv);
     58    if ((slit_trace_blu = FindSpline (argv[N])) == NULL) return (FALSE);
     59    remove_argument (N, &argc, argv);
     60  }
     61  if ((N = get_argument (argc, argv, "-psf-trace"))) {
     62    remove_argument (N, &argc, argv);
     63    if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE);
    3764    remove_argument (N, &argc, argv);
    3865  }
    … …  
    81108  }
    82109
     110  int VERBOSE = FALSE;
     111  if ((N = get_argument (argc, argv, "-v"))) {
     112    remove_argument (N, &argc, argv);
     113    VERBOSE = TRUE;
     114  }
     115
    83116  if (argc != 3) {
    84     gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-Nwave N] [-object vector] [-sky vector] [-backgnd vector] [-trace spline] [-profile vector]\n");
     117    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");
     118    gprint (GP_ERR, "  slit_trace_red : spline fit of slit central x pos vs y-coord for red chip\n");
     119    gprint (GP_ERR, "  slit_trace_blu : spline fit of slit central x pos vs y-coord for blu chip\n");
     120    gprint (GP_ERR, "  psf_trace      : spline fit of psf central x pos vs y-coord\n");
     121    gprint (GP_ERR, "  profile        : slit window profile (vector)\n");
     122    gprint (GP_ERR, "  object         : vector of object flux vs y-coord \n");
     123    gprint (GP_ERR, "  sky            : vector of local sky signal vs y-coord\n");
     124    gprint (GP_ERR, "  backgnd        : vector of extra-slit background flux vs y-coord\n");
     125    gprint (GP_ERR, "  PSF            : point-spread function vector (flux normalized, x-dir)\n");
     126    gprint (GP_ERR, "  LSF            : sline-spread function vector (flux normalized, y-dir)\n");
     127    gprint (GP_ERR, "  stilt          : slit tilt response : 2D kernel? \n");
     128    gprint (GP_ERR, "  redlimit       : pixel break between red and blu chips\n");
     129    gprint (GP_ERR, "  Nwave          : number of wavelength pixels in output image\n");
    85130    return FALSE;
    86131  }
    … …  
    123168  // if we are appying a trace offset spline, we need to generate an output window which
    124169  // is Nx + the full swing of the trace, then window back down
    125   if (trace) {
     170  if (slit_trace_red && slit_trace_blu) {
    126171    float dXmin = +1000;
    127172    float dXmax = -1000;
    128173    for (int iy = 0; iy < Ny; iy++) {
    129174      // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    130       float dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
     175      Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu;
     176      float dx = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
    131177      dXmin = MIN (dx, dXmin);
    132178      dXmax = MAX (dx, dXmax);
    … …  
    181227    if (psf) Npof -= Npsf;
    182228
     229    // integral and fractional pixel offsets of PSF (none, by default)
     230    int dxi = 0;
     231    float dxf = 0.0, dxr = 1.0;
     232   
     233    if (psf_trace) {
     234      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
     235      dxi = floor(dx);
     236      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
     237      dxr = 1 - dxf;
     238      Npof += dxi;
     239      // fprintf (stderr, "%d: %f : %d\n", iy, dx, Npof);
     240      if (VERBOSE && (iy % 200 == 0)) { gprint (GP_ERR, "%d : %f : %d : %f %f\n", iy, dx, dxi, dxf, dxr); }
     241    }
     242
     243    // if fractional offset is small, do not interpolate
     244    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
     245
    183246    // flux = obj * PSF + sky
    184247    for (int ix = 0; ix < Nx; ix++) {
    … …  
    191254        // only add in the flux if we are in range of the PSF
    192255        if ((n >= 0) && (n < NpsfFull)) {
    193           value += objVy * psfV[n];
     256          if (doInterp) {
     257            if ((n > 0) && (n < NpsfFull)) {
     258              if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
     259                value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
     260              }
     261            }
     262            if (n == 0) {
     263              if (isfinite(psfV[n])) {
     264                value += objVy*psfV[n]*dxr;
     265              }
     266            }
     267          } else {
     268            if (isfinite(psfV[n])) {
     269              value += objVy*psfV[n];
     270            }
     271          }
    194272        }
    195273      } else {
    … …  
    329407
    330408  // shift pixels in x based on the trace
    331   if (trace) {
     409  if (slit_trace_red && slit_trace_blu) {
    332410    ALLOCATE_PTR (outTraceBuffer, char, NxBase*Ny*sizeof(float));
    333411    float *outTrace = (float *) outTraceBuffer;
    … …  
    337415
    338416      // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    339       float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
     417      Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu;
     418      float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
    340419     
    341420      // extract the integer pixel offset and the fractional offset

  • trunk/Ohana/src/opihi/mana/deimos_mkslit.c

    r40659 r41159  
    77  // input parameters:
    88  //   slit profile response : vector of fractional flux vs x-coord
    9   //   trace                 : spline fit of slit central x pos vs y-coord
    10   //   flux                  : vector of input signal vs y-coord
     9  //   trace_ref,trace_blu   : spline fit of slit central x pos vs y-coord
     10  //   flux (optional)       : vector of input signal vs y-coord
     11  //   sky  (optional)       : vector of input signal vs y-coord
     12  //   fluxbins (optional)   : size of output image in y-direction
     13  //   NOTE: one of flux, sky, or fluxbins must be defined
     14
     15  // output : buff (an image with Nx defined by profile and Ny defined
     16  //                by flux, sky, or fluxbins)
    1117
    1218  int N;
    1319
    14   Vector *profile = NULL;
    15   Vector *flux    = NULL;
    16   Vector *sky     = NULL;
    17   Spline *trace   = NULL;
    18   Buffer *buff    = NULL;
     20  Vector *profile   = NULL;
     21  Vector *flux      = NULL;
     22  Vector *sky       = NULL;
     23  Spline *trace_red = NULL;
     24  Spline *trace_blu = NULL;
     25  Buffer *buff      = NULL;
    1926
    2027  // user-specified flux value (otherwise assumed to be 1)
    … …  
    4047  }
    4148
    42   if (argc != 4) {
    43     gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
     49  // either flux or fluxbins must be specified to define output size
     50  int profilebins = 0;
     51  if ((N = get_argument (argc, argv, "-profilebins"))) {
     52    remove_argument (N, &argc, argv);
     53    profilebins = atoi (argv[N]);
     54    remove_argument (N, &argc, argv);
     55  }
     56
     57  // for a red vs blu spline, we need to specify the split point
     58  // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
     59  // this more generic (need to define the ranges somewhere)
     60  int redlimit = 4096;
     61  if ((N = get_argument (argc, argv, "-redlimit"))) {
     62    remove_argument (N, &argc, argv);
     63    redlimit = atoi (argv[N]);
     64    remove_argument (N, &argc, argv);
     65  }
     66
     67  if (argc != 5) {
     68    gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace_red) (trace_blu) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
    4469    return FALSE;
    4570  }
    … …  
    5176
    5277  // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
    53   if ((profile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
    54   if ((trace   = FindSpline (argv[2])) == NULL) return (FALSE);
    55   if ((buff    = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE);
     78  if ((profile   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     79  if ((trace_red = FindSpline (argv[2])) == NULL) return (FALSE);
     80  if ((trace_blu = FindSpline (argv[3])) == NULL) return (FALSE);
     81  if ((buff      = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
    5682
    5783  // define the output window
    58   int Nx = 2*profile[0].Nelements + 1;
     84  // if no profile half-width is specified, use supplied profile
     85  int Nx = profilebins ? 2*profilebins + 1 : 2*profile[0].Nelements + 1;
    5986  int Ny = flux ? flux[0].Nelements : fluxbins;
    6087 
    … …  
    73100
    74101    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
     102    Spline *trace = (iy < redlimit) ? trace_red : trace_blu;
    75103    float dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
    76104   
    … …  
    90118
    91119      // set the sky in the entire window
    92       out[ix + iy*Nx] = sky ? sky->elements.Flt[iy] : 0.0;
     120      out[ix + iy*Nx] = NAN;
    93121
    94122      // equivalent coord in the profile:
    … …  
    98126
    99127      // a default value:
     128      float vsky = sky ? sky->elements.Flt[iy] : 0.0;
    100129      float vout = NAN;
    101130
    … …  
    121150      if (flux) vout *= flux->elements.Flt[iy];
    122151
    123       out[ix + iy*Nx] += vout;
     152      out[ix + iy*Nx] = vsky + vout;
    124153    }
    125154  }
Note: See TracChangeset for help on using the changeset viewer.