Changeset 41159 for trunk/Ohana/src/opihi/mana/deimos_mkmodel.c

Timestamp:

Nov 27, 2019, 11:04:51 AM (7 years ago)

Author:

eugene

Message:

improvments to deimos extraction tools

File:

• trunk/Ohana/src/opihi/mana/deimos_mkmodel.c (modified) (7 diffs)

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

@@ -86 +86 @@
 // generate a subimage of the full model for the range row - row + Ny
 // obj,sky,bck are subset vectors corresponding to the range row to row+Ny
-float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *trace, int Nx, int Ny, int row) {
-
-  float *raw = deimos_make_straight_image (obj, sky, psf, Nx, Ny, row);
+// the slit_trace_* and psf_trace refer to the full image starting at row
+float *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) {
+
+  float *raw = deimos_make_straight_image (obj, sky, psf, psf_trace, Nx, Ny, row);
   float *new = deimos_apply_tilt (raw, Nx, Ny);
   deimos_apply_profile (profile, new, Nx, Ny);
   deimos_add_background (bck, new, Nx, Ny, row);
-  float *out = deimos_apply_trace (trace, new, Nx, Ny, row);
+  float *out = deimos_apply_trace (slit_trace_red, slit_trace_blu, redlimit, new, Nx, Ny, row);
 
   free (raw);
@@ -104 +105 @@
 // *** generate a raw image with no slit tilt and no trace offsets.  apply the PSF to the
 // object flux, add in the sky flux
-float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row) {
+float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row) {
   OHANA_UNUSED_PARAM(row);
   
@@ -113 +114 @@
 
   // Xref is a global supplied by the user (defaults to Nx/2 if < 0)
-  int Xoff = (int)(Xref - Npsf);
+  int XoffRef = (int)(Xref - Npsf);
 
   // Nx : width of output window
@@ -123 +124 @@
     
     // we are generating the image for just the row range row to row+Ny
+    // note that obj & sky are subset vectors for just this range of pixels
     opihi_flt objVy = obj[iy]; // if we shift to using the row offset: obj[iy + row]
     opihi_flt skyVy = sky[iy]; // if we shift to using the row offset: sky[iy + row]
       
+    // integral and fractional pixel offsets of PSF (none, by default)
+    int dxi = 0;
+    float dxf = 0.0, dxr = 1.0;
+    int Xoff = XoffRef;
+    
+    if (psf_trace) {
+      // the psf_trace is referenced against the full image, so we need to add row to iy:
+      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy + row);
+      dxi = floor(dx);
+      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
+      dxr = 1 - dxf;
+      Xoff += dxi;
+      // fprintf (stderr, "%d: %f : %d\n", iy + row, dx, Xoff);
+    }
+
+    // if fractional offset is small, do not interpolate
+    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
+
     // flux = obj * PSF + sky
     for (int ix = 0; ix < Nx; ix++) {
@@ -136 +156 @@
       // only add in the flux if we are in range of the PSF
       if ((n >= 0) && (n < NpsfFull)) {
-        value += objVy * psfV[n];
+        if (doInterp) {
+          if ((n > 0) && (n < NpsfFull)) {
+            if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
+              value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
+            }
+          }
+          if (n == 0) {
+            if (isfinite(psfV[n])) {
+              value += objVy*psfV[n]*dxr;
+            }
+          }
+        } else {
+          if (isfinite(psfV[n])) {
+            value += objVy * psfV[n];
+          }
+        }
       }
       out[ix + iy*Nx] = value;
@@ -219 +254 @@
 
 // shift pixels in x based on the trace
-float *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row) {
+float *deimos_apply_trace (Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, float *input, int Nx, int Ny, int row) {
 
   ALLOCATE_PTR (output, float, Nx*Ny);
@@ -227 +262 @@
     // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
     // NOTE: spline is evaluated at full dispersion coordinate
-    float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy + row);
+    int iy_full = iy + row;
+    Spline *slit_trace = (iy_full < redlimit) ? slit_trace_red : slit_trace_blu;
+    float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy_full);
       
     // extract the integer pixel offset and the fractional offset