IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 40642 for trunk


Ignore:
Timestamp:
Mar 18, 2019, 5:02:59 PM (7 years ago)
Author:
eugene
Message:

cleanup opihi sort functions; add deimos fitobj function

Location:
trunk/Ohana/src/opihi
Files:
5 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.data/nnet_train.c

    r40371 r40642  
    33// these are local only
    44float nnet_sigmoid (float value);
    5 void  nnet_sortseq (float *X, int *Y, int N);
    65int   nnet_onelayer (Nnet *nnet, int L);
    76void  nnet_reset_Nabla (Nnet *nnet);
    … …  
    191190    // generate a random sequence : used to select random mini batches
    192191    for (int i = 0; i < Ntrial; i++) { seq[i] = i; rnd[i] = drand48(); }
    193     nnet_sortseq (rnd, seq, Ntrial);
     192    sort_float_index (rnd, seq, Ntrial);
    194193
    195194    int Npass = Ntrial / Nmini;
    … …  
    514513  return 1.0 / (1.0 + exp(-value));
    515514}
    516 
    517 void nnet_sortseq (float *X, int *Y, int N) {
    518 
    519 # define SWAPFUNC(A,B){ float ftmp; int itmp;   \
    520     ftmp = X[A]; X[A] = X[B]; X[B] = ftmp;      \
    521     itmp = Y[A]; Y[A] = Y[B]; Y[B] = itmp;      \
    522   }
    523 # define COMPARE(A,B)(X[A] < X[B])
    524 
    525   OHANA_SORT (N, COMPARE, SWAPFUNC);
    526 
    527 # undef SWAPFUNC
    528 # undef COMPARE
    529 
    530 }
    531 

  • trunk/Ohana/src/opihi/cmd.data/sort.c

    r40398 r40642  
    11# include "data.h"
    2 
    3 void fsortindex (opihi_flt *X, int *IDX, int N) {
    4 
    5 # define SWAPFUNC(A,B){ opihi_flt tmp; int itmp;        \
    6   tmp = X[A]; X[A] = X[B]; X[B] = tmp; \
    7   itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \
    8 }
    9 
    10 // # define COMPARE(A,B)(X[A] < X[B])
    11 # define COMPARE(A,B)((!isfinite(X[A]) && isfinite(X[B])) || (X[A] < X[B]))
    12 
    13   OHANA_SORT (N, COMPARE, SWAPFUNC);
    14 
    15 # undef SWAPFUNC
    16 # undef COMPARE
    17 
    18 }
    19 
    20 void isortindex (opihi_int *X, int *IDX, int N) {
    21 
    22 # define SWAPFUNC(A,B){ opihi_int tmp; int itmp;        \
    23   tmp = X[A]; X[A] = X[B]; X[B] = tmp; \
    24   itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \
    25 }
    26 
    27 # define COMPARE(A,B)(X[A] < X[B])
    28 // # define COMPARE(A,B)((!isfinite(X[A]) && isfinite(X[B])) || (X[A] < X[B]))
    29 
    30   OHANA_SORT (N, COMPARE, SWAPFUNC);
    31 
    32 # undef SWAPFUNC
    33 # undef COMPARE
    34 
    35 }
    362
    373// XXX add an option to NOT sort, but return an index instead?
    … …  
    7844  /* sort key & index */
    7945  if (vec[0][0].type == OPIHI_FLT) {
    80     fsortindex (vec[0][0].elements.Flt, index, Nval);
     46    sort_opihi_flt_index (vec[0][0].elements.Flt, index, Nval);
    8147  } else {
    82     isortindex (vec[0][0].elements.Int, index, Nval);
     48    sort_int_index (vec[0][0].elements.Int, index, Nval);
    8349  }
    8450

  • trunk/Ohana/src/opihi/dvo/gstar.c

    r40291 r40642  
    7171void printPhotcodeSequence (Average *average, SecFilt *secfilt, int entry, int type);
    7272
    73 void sort_index_float (float *X, int *IDX, int N) {
    74 
    75 # define SWAPFUNC(A,B){ float tmp; int itmp;            \
    76     tmp = X[A]; X[A] = X[B]; X[B] = tmp;                \
    77     itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp;      \
    78   }
    79 # define COMPARE(A,B)(X[A] < X[B])
    80 
    81   OHANA_SORT (N, COMPARE, SWAPFUNC);
    82 
    83 # undef SWAPFUNC
    84 # undef COMPARE
    85 
    86 }
    87 
    8873int gstar (int argc, char **argv) {
    8974 
    … …  
    712697        }
    713698
    714         sort_index_float (value, index, catalog.average[k].Nmeasure);
     699        sort_float_index (value, index, catalog.average[k].Nmeasure);
    715700
    716701        for (j = 0; j < catalog.average[k].Nmeasure; j++, m++) {

  • trunk/Ohana/src/opihi/include/data.h

    r40376 r40642  
    249249int GetKapaChannelFromString (char *string);
    250250
     251// in sort_funcs.c:
     252void sort_opihi_flt_index (opihi_flt *X, int *IDX, int N);
     253void sort_int_index (opihi_int *X, int *IDX, int N);
     254void sort_float_index (float *X, int *IDX, int N);
     255
    251256/*** Neural Network functions (nnet_*) ***/
    252257

  • trunk/Ohana/src/opihi/lib.data/Makefile

    r40319 r40642  
    2929$(SDIR)/nnet.$(ARCH).o                  \
    3030$(SDIR)/precess.$(ARCH).o               \
     31$(SDIR)/sort_funcs.$(ARCH).o            \
    3132$(SDIR)/starfuncs.$(ARCH).o             \
    3233$(SDIR)/hermitian.$(ARCH).o             \

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

    r40622 r40642  
    2424$(SRC)/deimos_mkslit.$(ARCH).o \
    2525$(SRC)/deimos_fitslit.$(ARCH).o \
     26$(SRC)/deimos_getobj.$(ARCH).o \
     27$(SRC)/deimos_fitobj.$(ARCH).o \
     28$(SRC)/deimos_mkalt.$(ARCH).o \
     29$(SRC)/deimos_mkmodel.$(ARCH).o \
    2630$(SRC)/findrowpeaks.$(ARCH).o
    2731

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

    r40622 r40642  
    44int deimos_fitslit (int argc, char **argv);
    55int deimos_mkobj (int argc, char **argv);
     6int deimos_mkalt (int argc, char **argv);
     7int deimos_getobj (int argc, char **argv);
     8int deimos_fitobj (int argc, char **argv);
    69
    710static Command deimos_commands[] = {
     11  {1, "fitobj", deimos_fitobj, "fit for object parameters"},
     12  {1, "getobj", deimos_getobj, "determine crude object parameters"},
    813  {1, "mkobj", deimos_mkobj, "make a full object image"},
     14  {1, "mkalt", deimos_mkalt, "make a full object image (uses deimos_make_object, for a test)"},
    915  {1, "mkslit", deimos_mkslit, "make a slit image"},
    1016  {1, "fitslit", deimos_fitslit, "fit slit image to observed slit flux"},

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

    r40640 r40642  
    11# include "data.h"
     2
     3// temporary storage used in local functions:
     4
     5float *tmpv = NULL;
     6float *tmpx = NULL;
     7
     8// pass in a pointer to the row (&buffV[iy*Nx])
     9float get_segment_median (float *buffV, float Xs, float Xe, float *Xref) {
     10
     11  // get the pixels in the left portion below profile:
     12  // XXX: worry about the effect of a fractional offset
     13  for (int ix = Xs + 0.5; ix < Xe + 0.5; ix++) {
     14    int i = ix - Xs;
     15    tmpv[i] = buffV[ix];
     16    tmpx[i] = ix;
     17  }
     18  int Nv = Xe - Xs;
     19
     20  // generate a statistic from these pixels: mean, median, etc?
     21  fsort (tmpv, Nv);
     22  int Nv2 = Nv/2; // Nv == 1, Nv2 = 0, Nv == 2, Nv2 = 1, Nv == 3, Nv2 = 1, Nv == 4, Nv == 2
     23  float V = (Nv % 2) ? tmpv[Nv2] : 0.5*(tmpv[Nv2-1] + tmpv[Nv2]);
     24   
     25  // calculate average coordinate (not flux-weighted...)
     26  // XXX : if we allow (and exclude) NAN-valued pixels, we can calculate this more explicitly
     27  *Xref = 0.5*(Xe + Xs);
     28
     29  return V;
     30}
    231
    332int deimos_getobj (int argc, char **argv) {
    … …  
    74103  } else { goto usage; }
    75104
    76   if (argc != 2) goto usage;
    77 
     105  if (argc != 3) goto usage;
     106
     107  // buffer is the full image, Xref is reference coordinate of the profile in the buffer
    78108  if ((buffer = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     109  float Xref = atof(argv[2]);
     110 
     111  // profile and PSF are defined with central reference pixel mapped to Xref (+ trace)
     112  int Nprofile = profile->Nelements;
     113  if (Nprofile % 2 == 0) {
     114    gprint (GP_ERR, "slit profile vector must have an odd number of pixels\n");
     115    return FALSE;
     116  }
     117  int Xprofile = Nprofile / 2;
     118
     119  int Npsf = psf->Nelements;
     120  if (Npsf % 2 == 0) {
     121    gprint (GP_ERR, "PSF vector must have an odd number of pixels\n");
     122    return FALSE;
     123  }
     124  int Xpsf = Npsf / 2;
     125
     126  // we will generate output vectors (object, sky, background) of length Ny
     127  int Nx = buffer[0].matrix.Naxis[0];
     128  int Ny = buffer[0].matrix.Naxis[1];
     129  ResetVector (object,  OPIHI_FLT, Ny);
     130  ResetVector (sky,     OPIHI_FLT, Ny);
     131  ResetVector (backgnd, OPIHI_FLT, Ny);
    79132 
    80133  // I need to generate a sky-left and sky-right vector in a region on the slit portion of
    … …  
    82135  // based on stilt.  Use these sky vectors to predict the sky under the PSF
    83136
    84   // Here are the code stages:
    85 
    86   // *********
    87 
    88   opihi_flt *objV = (opihi_flt *) (object ? object->elements.Flt : NULL);
    89   opihi_flt *skyV = (opihi_flt *) (sky    ? sky->elements.Flt : NULL);
    90   opihi_flt *psfV = (opihi_flt *) (psf    ? psf->elements.Flt : NULL);
    91   opihi_flt *lsfV = (opihi_flt *) (lsf    ? lsf->elements.Flt : NULL);
    92 
    93   // psf vector must have odd number of pixels:
    94   int Npsf = 0;
    95   int NpsfFull = 0;
    96   if (psf) {
    97     if (psf->Nelements % 2 == 0) {
    98       gprint (GP_ERR, "psf vector must have an odd number of pixels\n");
    99       return FALSE;
    100     }
    101     NpsfFull = psf->Nelements;
    102     Npsf = NpsfFull / 2;
    103   }
    104 
    105   // lsf vector must have odd number of pixels:
    106   int Nlsf = 0;
    107   int NlsfFull = 0;
    108   if (lsf) {
    109     if (lsf->Nelements % 2 == 0) {
    110       gprint (GP_ERR, "lsf vector must have an odd number of pixels\n");
    111       return FALSE;
    112     }
    113     NlsfFull = lsf->Nelements;
    114     Nlsf = NlsfFull / 2;
    115   }
    116 
    117   // *** apply the PSF to the object flux, add in the sky flux
    118 
    119   ALLOCATE_PTR (s1, float, Nx*Ny);
     137  // *** scan the profile to determine the background and sky cut points ***
     138
     139  // values below in the profile coordinate system need to be shifted by (Xref - Nprofile/2)
     140  // to match the image coordinate system
     141
     142  int B_s_l = 0, B_e_l = -1; // left background window
     143  int B_s_r = -1, B_e_r = Nprofile - 1; // right background window
     144  int S_left = -1, S_right = -1; // top portion of profile
     145  for (int i = 0; i < Nprofile; i++) {
     146    if ((B_e_l   <  0) && (profile->elements.Flt[i] > 0.10)) B_e_l = i - 1;
     147    if ((B_e_l   >= 0) && (S_left  < 0) && (profile->elements.Flt[i] > 0.95)) S_left = i;
     148    if ((S_left  >= 0) && (S_right < 0) && (profile->elements.Flt[i] < 0.95)) S_right = i - 1;
     149    if ((S_right >= 0) && (profile->elements.Flt[i] < 0.10)) { B_s_r = i; break; }
     150  }
     151  if (B_e_l == B_s_l) B_e_l ++;
     152  if (B_e_r == B_s_r) B_s_r --;
     153
     154  // S_left to S_right defines are window at the top, choose 0.05 - 0.15, 0.85 - 0.95
     155  int N_top = S_right - S_left + 1;
     156
     157  // if N_top < 20, these jumps might be too small
     158  // if N_top < 20, these jumps might be too small
     159  int S_s_l = S_left + 0.05*N_top;
     160  int S_e_l = S_left + 0.15*N_top;
     161  if (S_e_l == S_s_l) S_e_l ++;
     162  int S_s_r = S_left + 0.85*N_top;
     163  int S_e_r = S_left + 0.95*N_top;
     164  if (S_e_r == S_s_r) S_s_r --;
     165
     166  gprint (GP_ERR, "windows: %.1f - %.1f, %.1f - %.1f, %.1f - %.1f, %.1f - %.1f\n",
     167          B_s_l + Xref - Xprofile, B_e_l + Xref - Xprofile,
     168          S_s_l + Xref - Xprofile, S_e_l + Xref - Xprofile,
     169          S_s_r + Xref - Xprofile, S_e_r + Xref - Xprofile,
     170          B_s_r + Xref - Xprofile, B_e_r + Xref - Xprofile);
     171
     172  // OPEN QUESTIONS:
     173  // NOTE: if the range of pixels in the windows is too large, then the tilt will wash out the signal.
     174  // this pass is a guess, so perhaps that does not matter, but there is clearly a trade-off here.
     175
     176  // generate two buffers to store the flux and coordinate for the background, sky, object pixels:
     177  ALLOCATE (tmpv, float, Nx);
     178  ALLOCATE (tmpx, float, Nx);
     179
     180  float *buffV = (float *) buffer[0].matrix.buffer;
     181
     182  // use these to store the flux vectors
     183  ALLOCATE_PTR (bckF_lf, float, Ny);
     184  ALLOCATE_PTR (bckF_rt, float, Ny);
     185  ALLOCATE_PTR (skyF_lf, float, Ny);
     186  ALLOCATE_PTR (skyF_rt, float, Ny);
     187
     188  // use these to store the x-coord vectors
     189  ALLOCATE_PTR (bckX_lf, float, Ny);
     190  ALLOCATE_PTR (bckX_rt, float, Ny);
     191  ALLOCATE_PTR (skyX_lf, float, Ny);
     192  ALLOCATE_PTR (skyX_rt, float, Ny);
     193
     194  float dxMax = 0; // maximum distance between sky regions
    120195
    121196  // loop over the y positions
    122197  for (int iy = 0; iy < Ny; iy++) {
    123198
    124     opihi_flt objVy = objV ? objV[iy] : 1.0; // if object is not supplied, flux defaults to 1.0
    125     opihi_flt skyVy = skyV ? skyV[iy] : 0.0; // if sky is not supplied, flux defaults to 0.0
     199    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
     200    float dXtrace = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
     201
     202    // Xref is the reference coordinate of the profile on the image
     203    // Xprofile is the center of the profile. 
     204    // dXtrace follows the x-dir variations of the slit profile
     205    float Xtrace = Xref - Xprofile + dXtrace;
     206
     207    bckF_lf[iy] = get_segment_median(&buffV[iy*Nx], B_s_l + Xtrace, B_e_l + Xtrace, &bckX_lf[iy]);
     208    bckF_rt[iy] = get_segment_median(&buffV[iy*Nx], B_s_r + Xtrace, B_e_r + Xtrace, &bckX_rt[iy]);
     209    skyF_lf[iy] = get_segment_median(&buffV[iy*Nx], S_s_l + Xtrace, S_e_l + Xtrace, &skyX_lf[iy]);
     210    skyF_rt[iy] = get_segment_median(&buffV[iy*Nx], S_s_r + Xtrace, S_e_r + Xtrace, &skyX_rt[iy]);
     211    dxMax = MAX (dxMax, skyX_rt[iy] - skyX_lf[iy]); // this is a constant value unless we ignore some of the pixels
     212    // assert skyX_rt > skyX_lf ?
     213  }
     214
     215  // extract sky-subtracted values for the psf fit, using the two sky vectors to predict the sky
     216  float dyMax = dxMax*sin(stilt*RAD_DEG);  // dxMax is based on the profile width, calculated above
     217
     218  // find nominal dy_lf, dy_rt offsets for the sky vectors. 
     219  int dy_lf = -dyMax / 2;
     220  int dy_rt = +dyMax / 2;
     221
     222  // we are fitting this PSF to the data:
     223  opihi_flt *psfV = psf->elements.Flt;
     224
     225  // init the regions which are outside the valid region
     226  for (int iy = 0; iy < fabs(dyMax); iy++) {
     227    object->elements.Flt[iy] = 0;
     228    sky->elements.Flt[iy] = 0;
     229    backgnd->elements.Flt[iy] = 0;
     230  }   
     231  for (int iy = Ny - fabs(dyMax); iy < Ny; iy++) {
     232    object->elements.Flt[iy] = 0;
     233    sky->elements.Flt[iy] = 0;
     234    backgnd->elements.Flt[iy] = 0;
     235  }   
     236
     237
     238  // Ys, Ye based on the dyMax value above (skip the first few lines to avoid running out of bounds on sky)
     239  for (int iy = fabs(dyMax); iy < Ny - fabs(dyMax); iy++) {
     240
     241    // X_lf & X_rt are calculated in image x-coords
     242    float S_lf = skyF_lf[iy + dy_lf];
     243    float S_rt = skyF_rt[iy + dy_rt];
     244    float X_lf = skyX_lf[iy + dy_lf];
     245    float X_rt = skyX_rt[iy + dy_rt];
     246   
     247    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
     248    float dXtrace = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
    126249     
    127     int Npof = floor(Nx/2);
    128     if (psf) Npof -= Npsf;
    129 
    130     // flux = obj * PSF + sky
    131     for (int ix = 0; ix < Nx; ix++) {
    132      
    133       opihi_flt value = skyVy;
    134 
    135       if (psfV) {
    136         int n = ix - Npof;
    137         // n is the pixel in the PSF corresponding to the ix pixel in the output image
    138         // only add in the flux if we are in range of the PSF
    139         if ((n >= 0) && (n < NpsfFull)) {
    140           value += objVy * psfV[n];
    141         }
    142       } else {
    143         if (ix == Nx / 2) value += objVy;
    144       }   
    145       s1[ix + iy*Nx] = value;
     250    // Xref is the reference coordinate of the profile on the image
     251    // Xpsf is the center of the PSF. 
     252    // dXtrace follows the x-dir variations of the slit profile (or PSF?)
     253    float Xtrace = Xref - Xpsf + dXtrace;
     254   
     255    // now I need to find the limits of the PSF in the buffer at this iy location
     256    for (int ipsf = 0; ipsf < Npsf; ipsf ++) {
     257
     258      int ix = ipsf + Xtrace + 0.5; // image coordinate (nearest pixel center)
     259      float So = S_lf + (ix - X_lf)*(S_rt - S_lf)/(X_rt - X_lf); // interpolated sky under PSF peak
     260
     261      tmpv[ipsf] = buffV[ix + iy*Nx] - So;
     262      tmpx[ipsf] = ipsf;
    146263    }
    147   }
    148 
    149   // *** apply the slit tilt & lsf
    150 
    151   // first, generate the interpolation kernels.  For a slit tilt of theta, there is a
    152   // displacement in the y-direction of dy = dx sin(theta) where dx is the x-coord
    153   // relative to the slit window center (Nx / 2).  we thus need an image of size Nx,
    154   // Nx*sin(theta)
    155 
    156   float sin_stilt = sin(stilt*RAD_DEG);
    157   int Nyk = ceil(Nx * fabs(sin_stilt));
    158   if (Nyk < 3) Nyk = 3; // minimum of 3 pixels (or kernel assumption below fails)
    159   Nyk += NlsfFull;
    160   if (Nyk % 2 == 0) Nyk ++; // force Nyk to be odd
    161   int Nyk2 = Nyk/2; // Nyk2 is the center pixel of the interpolations
    162   ALLOCATE_PTR (kernel, float, Nx*Nyk);
    163   int Nkpix = Nx*Nyk;
    164 
    165   int Nx2 = floor(Nx/2);
    166   for (int ix = 0; ix < Nx; ix++) {
    167     for (int iy = 0; iy < Nyk; iy++) { kernel[ix + iy*Nx] = 0.0; }
    168 
    169     // displacement in y-dir due to slit tilt:
    170     float dy = (ix - Nx2) * sin_stilt;
    171     int dyi = floor(dy);
    172     float dyf = dy - dyi;
    173 
    174     if (lsf) {
    175       int Sy = Nyk2 - Nlsf + dyi;
    176       int doInterp = (fabs(dyf) < 1e-5) ? FALSE : TRUE;
    177       for (int iy = 0; iy < Nyk; iy++) {
    178         int py = iy - Sy;
    179         if (py < 0) continue;
    180         if (py >= NlsfFull) continue;
    181 
    182         float vout = NAN;
    183         if (doInterp) {
    184           if ((py > 0) && (py < NlsfFull - 1)) {
    185             vout = lsfV[py]*(1 - dyf) + lsfV[py-1]*dyf;
    186           }
    187           if (py == 0) {
    188             vout = lsfV[py]*dyf;
    189           }
    190           if (py == NlsfFull - 1) {
    191             vout = lsfV[py]*(1.0 - dyf);
    192           }
    193         } else {
    194           vout = lsfV[py];
    195         }
    196         int Nkpix_i = ix + iy*Nx;
    197         myAssert (Nkpix_i < Nkpix, "oops");
    198         kernel[Nkpix_i] = vout;
    199       }
    200     } else {
    201       int Nkpix_i;
    202       Nkpix_i = ix + (dyi + Nyk2 + 0)*Nx;
    203       myAssert (Nkpix_i < Nkpix, "oops");
    204       kernel[Nkpix_i] = 1.0 - dyf;
    205 
    206       Nkpix_i = ix + (dyi + Nyk2 + 1)*Nx;
    207       myAssert (Nkpix_i < Nkpix, "oops");
    208       kernel[Nkpix_i] = dyf;
    209       fprintf (stderr, "%d, %d: %f %f\n", ix, (dyi + Nyk2), kernel[ix + (dyi + Nyk2 + 0)*Nx], kernel[ix + (dyi + Nyk2 + 1)*Nx]);
     264
     265    // I now have a pair of vectors (tmpx, tmpv) which are aligned to the PSF at this location
     266    // now fit the PSF to tmpv by calculating a normalization
     267    float S1 = 0.0, S2 = 0.0;
     268    for (int ipsf = 0; ipsf < Npsf; ipsf++) {
     269      S1 += tmpv[ipsf]*psfV[ipsf]; // if we have an errorbar, include / SQ(sigv[i])
     270      S2 += psfV[ipsf]*psfV[ipsf]; // if we have an errorbar, include / SQ(sigv[i])
    210271    }
    211   }
    212 
    213   float *out = (float *) output[0].matrix.buffer;
    214 
    215   // next, apply the interpolation kernel to the image
    216 
    217   // loop over the y positions
    218   for (int iy = 0; iy < Ny; iy++) {
    219 
    220     // use the tilt kernel to interpolate to this x coordinate
    221     for (int ix = 0; ix < Nx; ix++) {
    222 
    223       // if the LSF is defined, we cannot bypass the convolution
    224       if (!lsf && (fabs(stilt) < 1e-5)) {
    225         out[ix + iy*Nx] = s1[ix + iy*Nx];
    226       } else {
    227         opihi_flt g = 0.0, s = 0.0;
    228         for (int n = -Nyk2; n <= Nyk2; n++) {
    229           if (iy + n < 0) continue; // bottom edge of full image
    230           if (iy + n >= Ny) continue; // top edge of full image
    231           int Nkpix_i = ix + (n + Nyk2)*Nx;
    232           myAssert (Nkpix_i < Nkpix, "oops");
    233           s += kernel[Nkpix_i]*s1[ix + (iy + n)*Nx];
    234           g += kernel[Nkpix_i];
    235         }
    236         out[ix + iy*Nx] = s / g;
    237       }
    238     }
    239   }
    240   free (kernel);
    241   free (s1);
    242  
    243   // next, apply the profile (output = input * profile)
    244   opihi_flt *profileV = (opihi_flt *) (profile ? profile->elements.Flt : NULL);
    245   if (profile) {
    246     int Nprof = profile->Nelements;
    247     int Nprof2 = Nprof / 2;
    248 
    249     // loop over the y positions
    250     for (int iy = 0; iy < Ny; iy++) {
    251      
    252       int Sx = (int)(Nx2 - Nprof2);
    253      
    254       for (int ix = 0; ix < Nx; ix++) {
    255        
    256         // equivalent coord in the profile:
    257         int px = ix - Sx;
    258         if ((px >= 0) && (px < Nprof)) {
    259           out[ix + iy*Nx] *= profileV[px];
    260         } else {
    261           out[ix + iy*Nx] = 0.0;
    262         }
    263       }
    264     }
    265   }
    266 
    267   // add in the background
    268   opihi_flt *backgndV = (opihi_flt *) (backgnd ? backgnd->elements.Flt : NULL);
    269   if (backgnd) {
    270     for (int iy = 0; iy < Ny; iy++) {
    271       for (int ix = 0; ix < Nx; ix++) {
    272         out[ix + iy*Nx] += backgndV[iy];
    273       }
    274     }
    275   }
    276 
    277   // shift pixels in x based on the trace
    278   if (trace) {
    279     ALLOCATE_PTR (outTraceBuffer, char, NxBase*Ny*sizeof(float));
    280     float *outTrace = (float *) outTraceBuffer;
    281     int NxBase2 = NxBase / 2;
    282 
    283     for (int iy = 0; iy < Ny; iy++) {
    284 
    285       // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    286       float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
    287      
    288       // extract the integer pixel offset and the fractional offset
    289       int dxi = floor(dx); // -1.7 -> -2, -0.5 -> -1, +0.5 -> 0, +1.7 -> 1
    290       float dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
    291       float dxr = 1 - dxf;
    292 
    293       // if fractional offset is small, do not interpolate
    294       int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
    295      
    296       // How do I handle this?  define a temporary window which is a fraction of the full output window?
    297       int Sx = Nx2 - NxBase2 + dxi;
    298  
    299       // save the original output row
    300       // for (int ix = 0; ix < Nx; ix++) { tmprow[ix] = out[ix + iy*Nx]; }
    301 
    302       for (int ix = 0; ix < NxBase; ix++) {
    303 
    304         // equivalent coord in temporary buffer (Nx * Ny):
    305         int px = ix + Sx;
    306         if (px < 0) {
    307           outTrace[ix + iy*NxBase] = out[iy*Nx];
    308           continue;
    309         }
    310         if (px >= Nx) {
    311           outTrace[ix + iy*NxBase] = out[(Nx - 1) + iy*Nx];
    312           continue;
    313         }
    314        
    315         // a default value:
    316         float vout = NAN;
    317        
    318         int Npix = px + iy*Nx;
    319 
    320         if (doInterp) {
    321           if ((px > 0) && (px < Nx - 1)) {
    322             vout = out[Npix]*dxr + out[Npix + 1]*dxf;
    323           }
    324           if (px == 0) {
    325             vout = out[Npix];
    326           }
    327           if (px == Nx - 1) {
    328             vout = out[Npix];
    329           }
    330         } else {
    331           vout = out[Npix];
    332         }
    333         outTrace[ix + iy*NxBase] = vout;
    334       }
    335     }
    336 
    337     ResetBuffer (output, NxBase, Ny, -32, 0.0, 1.0);
    338     free (output[0].matrix.buffer);
    339     output[0].matrix.buffer = outTraceBuffer;
    340   }
     272
     273    int Xo = Xref - dXtrace; // center of PSF in image x-coordinates
     274    float So = S_lf + (Xo - X_lf)*(S_rt - S_lf)/(X_rt - X_lf); // interpolated sky under PSF peak
     275    object->elements.Flt[iy] = S1 / S2;
     276    sky->elements.Flt[iy] = So;
     277    backgnd->elements.Flt[iy] = 0.5*(bckF_lf[iy] + bckF_rt[iy]);
     278  }
     279
     280  free (tmpv);
     281  free (tmpx);
     282
     283  free (bckF_lf);
     284  free (bckF_rt);
     285  free (skyF_lf);
     286  free (skyF_rt);
     287
     288  free (bckX_lf);
     289  free (bckX_rt);
     290  free (skyX_lf);
     291  free (skyX_rt);
    341292
    342293  return TRUE;
    343294
    344295 usage:
    345   gprint (GP_ERR, "USAGE: deimos getobj (buffer) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n");
     296  gprint (GP_ERR, "USAGE: deimos getobj (buffer) (Xref) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n");
    346297  // can I make default values for trace (0.0), profile (?), psf (fraction of profile?), stilt (0.0)
    347298  return FALSE;
Note: See TracChangeset for help on using the changeset viewer.