Changeset 40642 for trunk/Ohana/src/opihi/mana/deimos_getobj.c

Timestamp:

Mar 18, 2019, 5:02:59 PM (7 years ago)

Author:

eugene

Message:

cleanup opihi sort functions; add deimos fitobj function

File:

• trunk/Ohana/src/opihi/mana/deimos_getobj.c (modified) (3 diffs)

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

@@ -1 +1 @@
 # 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 +103 @@
   } 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 +135 @@
   // 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;