IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 31, 2019, 3:08:13 PM (7 years ago)
Author:
eugene
Message:

deimos updates to add psf trace

Location:
branches/eam_branches/ohana.20190329/src/opihi
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ohana.20190329/src/opihi/cmd.data/interpolate_presort.c

    r40007 r40663  
    77  double x0, dx, dy, y0;
    88  Vector *xout, *yout, *xin, *yin;
     9
     10  int N;
     11  int FillEnds = FALSE;
     12  if ((N = get_argument (argc, argv, "-fill-ends"))) {
     13    FillEnds = TRUE;
     14    remove_argument (N, &argc, argv);
     15  }
    916
    1017  /** check basic syntax **/
    … …  
    3845  /* in and out vectors are sorted */
    3946  j = 0;
     47  i = 0;
     48
     49  // fill in the rest with NANs
     50  while (xout[0].elements.Flt[i] < xin[0].elements.Flt[0]) {
     51    yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[0] : NAN;
     52    i++;
     53  }
    4054
    4155  // every output pixel should get a value unless
    4256  // we are below the lowest in or above the highest in
    43   for (i = 0, j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) {
     57  for (j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) {
    4458
    4559    yout[0].elements.Flt[i] = NAN;
    … …  
    6680
    6781  // fill in the rest with NANs
     82  int NinLast = xin[0].Nelements - 1;
    6883  while (i < yout[0].Nelements) {
    69     yout[0].elements.Flt[i] = NAN;
     84    yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[NinLast] : NAN;
    7085    i++;
    7186  }
    7287
    7388  return (TRUE);
    74    
    7589}

  • branches/eam_branches/ohana.20190329/src/opihi/include/deimos.h

    r40642 r40663  
    77} DeimosResult;
    88
    9 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);
     9float     *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace, Spline *psf_trace, int Nx, int Ny, int row);
    1010void       deimos_make_kernel (float stilt, int Nx);
    1111void       deimos_free_kernel ();
    … …  
    1313
    1414// internal to make_model:
    15 float     *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row);
     15float     *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row);
    1616float     *deimos_apply_tilt (float *input, int Nx, int Ny);
    1717void       deimos_apply_profile (Vector *profile, float *out, int Nx, int Ny);
    1818void       deimos_add_background (opihi_flt *backgnd, float *out, int Nx, int Ny, int row);
    19 float     *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row);
     19float     *deimos_apply_trace (Spline *slit_trace, float *input, int Nx, int Ny, int row);
    2020

  • branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile

    r40655 r40663  
    2525$(SRC)/deimos_fitslit.$(ARCH).o \
    2626$(SRC)/deimos_getobj.$(ARCH).o \
     27$(SRC)/deimos_getalt.$(ARCH).o \
    2728$(SRC)/deimos_fitobj.$(ARCH).o \
    2829$(SRC)/deimos_fitalt.$(ARCH).o \

  • branches/eam_branches/ohana.20190329/src/opihi/mana/deimos.c

    r40655 r40663  
    66int deimos_mkalt (int argc, char **argv);
    77int deimos_getobj (int argc, char **argv);
     8int deimos_getalt (int argc, char **argv);
    89int deimos_fitobj (int argc, char **argv);
    910int deimos_fitalt (int argc, char **argv);
    … …  
    1415  {1, "fitalt", deimos_fitalt, "fit for object parameters using LMM"},
    1516  {1, "getobj", deimos_getobj, "determine crude object parameters"},
     17  {1, "getalt", deimos_getalt, "determine crude object parameters"},
    1618  {1, "mkobj", deimos_mkobj, "make a full object image"},
    1719  {1, "mkalt", deimos_mkalt, "make a full object image (uses deimos_make_object, for a test)"},

  • branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitalt.c

    r40655 r40663  
    2727  // input parameters:
    2828  // * buffer      : 2D image of slit (full frame or cutout?)
    29   // * trace       : spline fit of slit central x pos vs y-coord
     29  // * slit_trace  : spline fit of slit central x pos vs y-coord
     30  // * psf_trace   : spline fit of slit central x pos vs y-coord
    3031  // * profile     : slit window profile (vector)
    3132  // * PSF         : point-spread function vector (flux normalized, x-dir)
    … …  
    3940  int N;
    4041
    41   Spline *trace   = NULL;
     42  Spline *slit_trace = NULL;
     43  Spline *psf_trace  = NULL;
    4244  Vector *profile = NULL;
    4345
    … …  
    106108
    107109  // Input parameters:
    108   if ((N = get_argument (argc, argv, "-trace"))) {
    109     remove_argument (N, &argc, argv);
    110     if ((trace = FindSpline (argv[N])) == NULL) return (FALSE);
     110  if ((N = get_argument (argc, argv, "-slit-trace"))) {
     111    remove_argument (N, &argc, argv);
     112    if ((slit_trace = FindSpline (argv[N])) == NULL) return (FALSE);
     113    remove_argument (N, &argc, argv);
     114  } else { goto usage; }
     115  if ((N = get_argument (argc, argv, "-psf-trace"))) {
     116    remove_argument (N, &argc, argv);
     117    if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE);
    111118    remove_argument (N, &argc, argv);
    112119  } else { goto usage; }
    … …  
    244251      // current vector value:
    245252      // for deimos_make_model, row is the starting point of the subimage
    246       float *model_ref = deimos_make_model (objCurr, skyCurr, bckCurr, psf, profile, trace, Nx, Nrow, row_0);
     253      float *model_ref = deimos_make_model (objCurr, skyCurr, bckCurr, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0);
    247254
    248255      // in the test functions below, dRow is the central pixel of the subset, e.g., obj[row] => objCurr[dRow]
    … …  
    252259      // delta obj (save dObj)
    253260      opihi_flt *objtest = deimos_make_test (objCurr, objNoise[row], dRow, Nrow, &dObj);
    254       float *model_obj = deimos_make_model (objtest, skyCurr, bckCurr, psf, profile, trace, Nx, Nrow, row_0);
     261      float *model_obj = deimos_make_model (objtest, skyCurr, bckCurr, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0);
    255262
    256263      // delta sky (save dSky)
    257264      opihi_flt *skytest = deimos_make_test (skyCurr, skyNoise[row], dRow, Nrow, &dSky);
    258       float *model_sky = deimos_make_model (objCurr, skytest, bckCurr, psf, profile, trace, Nx, Nrow, row_0);
     265      float *model_sky = deimos_make_model (objCurr, skytest, bckCurr, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0);
    259266
    260267      // delta bck (save dBck)
    261268      opihi_flt *bcktest = deimos_make_test (bckCurr, bckNoise[row], dRow, Nrow, &dBck);
    262       float *model_bck = deimos_make_model (objCurr, skyCurr, bcktest, psf, profile, trace, Nx, Nrow, row_0);
     269      float *model_bck = deimos_make_model (objCurr, skyCurr, bcktest, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0);
    263270
    264271      // dObj, dSky, dBck carry the delta used derivatives in and the delta applied to the values out
    … …  
    294301    }
    295302
    296     float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, trace, Nx, Ny, 0);
     303    float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, slit_trace, psf_trace, Nx, Ny, 0);
    297304
    298305    int Npts = 0;
    … …  
    305312
    306313 usage:
    307   gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n");
     314  gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -slit-trace spline -psf-trace spline -profile vector -psf vector -stilt (angle)\n");
    308315  return FALSE;
    309316}

  • branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_getobj.c

    r40659 r40663  
    33// temporary storage used in local functions:
    44
    5 float *tmpv = NULL;
    6 float *tmpx = NULL;
     5static float *tmpv = NULL;
     6static float *tmpx = NULL;
    77
    88// pass in a pointer to the row (&buffV[iy*Nx])
    … …  
    3939  // input parameters:
    4040  // * buffer      : 2D image of slit (full frame or cutout?)
    41   // * trace       : spline fit of slit central x pos vs y-coord
     41  // * slit_trace  : spline fit of slit central x pos vs y-coord
     42  // * psf_trace   : spline fit of psf central x pos vs y-coord
    4243  // * profile     : slit window profile (vector)
    4344  // * PSF         : point-spread function vector (flux normalized, x-dir)
    … …  
    5152  int N;
    5253
    53   Spline *trace   = NULL;
     54  Spline *slit_trace = NULL;
     55  Spline *psf_trace  = NULL;
    5456  Vector *profile = NULL;
    5557
    … …  
    7173
    7274  // Input parameters:
    73   if ((N = get_argument (argc, argv, "-trace"))) {
    74     remove_argument (N, &argc, argv);
    75     if ((trace = FindSpline (argv[N])) == NULL) return (FALSE);
     75  if ((N = get_argument (argc, argv, "-slit-trace"))) {
     76    remove_argument (N, &argc, argv);
     77    if ((slit_trace = FindSpline (argv[N])) == NULL) return (FALSE);
     78    remove_argument (N, &argc, argv);
     79  } else { goto usage; }
     80  if ((N = get_argument (argc, argv, "-psf-trace"))) {
     81    remove_argument (N, &argc, argv);
     82    if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE);
    7683    remove_argument (N, &argc, argv);
    7784  } else { goto usage; }
    … …  
    218225
    219226  // loop over the y positions
     227  // examine the slit window profile to measure background and slit sky
    220228  for (int iy = 0; iy < Ny; iy++) {
    221229
    222230    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    223     float dXtrace = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
     231    float dXtrace = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
    224232
    225233    // Xref is the reference coordinate of the profile on the image
    … …  
    284292   
    285293    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    286     float dXtrace = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
     294    float dXtraceSLIT = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
     295    float dXtracePSF  = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
    287296     
    288297    // Xref is the reference coordinate of the profile on the image
    289298    // Xpsf is the center of the PSF. 
    290299    // dXtrace follows the x-dir variations of the slit profile (or PSF?)
    291     float Xtrace = Xref - Xpsf + dXtrace;
    292    
     300    float XtracePSF = Xref - Xpsf + dXtraceSLIT + dXtracePSF;
     301    float XtraceSLIT = Xref - Xprofile + dXtraceSLIT;
     302
    293303    // include the PSF trace and restrict the fit below to portions where the profile is > 0.95
    294304
    295305    // now I need to find the limits of the PSF in the buffer at this iy location
     306    int Ntmp = 0;
    296307    for (int ipsf = 0; ipsf < Npsf; ipsf ++) {
    297308
    298       int ix = ipsf + Xtrace + 0.5; // image coordinate (nearest pixel center)
     309      int ix = ipsf + XtracePSF + 0.5; // image coordinate (nearest pixel center)
     310      int islit = ix - XtraceSLIT - 0.5; // image coordinate (nearest pixel center)
     311
    299312      float So = S_lf + (ix - X_lf)*(S_rt - S_lf)/(X_rt - X_lf); // interpolated sky under PSF peak
    300313
    301       tmpv[ipsf] = buffV[ix + iy*Nx] - So;
    302       tmpx[ipsf] = ipsf;
     314      if (islit < 0) continue;
     315      if (islit >= profile->Nelements) continue;
     316
     317      if (iy == 6000) {
     318        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]);
     319      }
     320
     321      if (profile->elements.Flt[islit] < 0.95) continue;
     322
     323      tmpv[Ntmp] = buffV[ix + iy*Nx] - So;
     324      tmpx[Ntmp] = psfV[ipsf];
     325      Ntmp ++;
    303326    }
    304327
    … …  
    306329    // now fit the PSF to tmpv by calculating a normalization
    307330    float S1 = 0.0, S2 = 0.0;
    308     for (int ipsf = 0; ipsf < Npsf; ipsf++) {
     331    for (int ipsf = 0; ipsf < Ntmp; ipsf++) {
    309332      if (iy == 6000) {
    310         fprintf (stderr, "%d : %f %f : %f\n", ipsf, tmpv[ipsf], psfV[ipsf], tmpv[ipsf]*psfV[ipsf]);
     333        fprintf (stderr, "%d : %f %f : %f\n", ipsf, tmpv[ipsf], tmpx[ipsf], tmpv[ipsf]*tmpx[ipsf]);
    311334      }
    312       S1 += tmpv[ipsf]*psfV[ipsf]; // if we have an errorbar, include / SQ(sigv[i])
    313       S2 += psfV[ipsf]*psfV[ipsf]; // if we have an errorbar, include / SQ(sigv[i])
     335      S1 += tmpv[ipsf]*tmpx[ipsf]; // if we have an errorbar, include / SQ(sigv[i])
     336      S2 += tmpx[ipsf]*tmpx[ipsf]; // if we have an errorbar, include / SQ(sigv[i])
    314337    }
    315338
    316     int Xo = Xref - dXtrace; // center of PSF in image x-coordinates
     339    int Xo = Xref - dXtraceSLIT - dXtracePSF; // center of PSF in image x-coordinates
    317340    float So = S_lf + (Xo - X_lf)*(S_rt - S_lf)/(X_rt - X_lf); // interpolated sky under PSF peak
    318341    object->elements.Flt[iy] = S1 / S2;

  • branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkmodel.c

    r40659 r40663  
    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);
     88float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace, Spline *psf_trace, int Nx, int Ny, int row) {
     89
     90  float *raw = deimos_make_straight_image (obj, sky, psf, psf_trace, Nx, Ny, row);
    9191  float *new = deimos_apply_tilt (raw, Nx, Ny);
    9292  deimos_apply_profile (profile, new, Nx, Ny);
    9393  deimos_add_background (bck, new, Nx, Ny, row);
    94   float *out = deimos_apply_trace (trace, new, Nx, Ny, row);
     94  float *out = deimos_apply_trace (slit_trace, new, Nx, Ny, row);
    9595
    9696  free (raw);
    … …  
    104104// *** generate a raw image with no slit tilt and no trace offsets.  apply the PSF to the
    105105// 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) {
     106float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row) {
    107107  OHANA_UNUSED_PARAM(row);
    108108 
    … …  
    126126    opihi_flt skyVy = sky[iy]; // if we shift to using the row offset: sky[iy + row]
    127127     
     128    // integral and fractional pixel offsets of PSF (none, by default)
     129    int dxi = 0;
     130    float dxf = 0.0, dxr = 1.0;
     131   
     132    if (psf_trace) {
     133      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
     134      dxi = floor(dx);
     135      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
     136      dxr = 1 - dxf;
     137      Xref += dxi;
     138    }
     139
     140    // if fractional offset is small, do not interpolate
     141    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
     142
    128143    // flux = obj * PSF + sky
    129144    for (int ix = 0; ix < Nx; ix++) {
    … …  
    136151      // only add in the flux if we are in range of the PSF
    137152      if ((n >= 0) && (n < NpsfFull)) {
    138         value += objVy * psfV[n];
     153        if (doInterp) {
     154          if ((n > 0) && (n < NpsfFull)) {
     155            value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
     156          }
     157          if (n == 0) {
     158            value += objVy*psfV[n]*dxr;
     159          }
     160        } else {
     161          value += objVy * psfV[n];
     162        }
    139163      }
    140164      out[ix + iy*Nx] = value;
    … …  
    219243
    220244// shift pixels in x based on the trace
    221 float *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row) {
     245float *deimos_apply_trace (Spline *slit_trace, float *input, int Nx, int Ny, int row) {
    222246
    223247  ALLOCATE_PTR (output, float, Nx*Ny);
    … …  
    227251    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    228252    // NOTE: spline is evaluated at full dispersion coordinate
    229     float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy + row);
     253    float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy + row);
    230254     
    231255    // extract the integer pixel offset and the fractional offset
Note: See TracChangeset for help on using the changeset viewer.