Changeset 40624

Timestamp:

Feb 23, 2019, 10:23:16 AM (7 years ago)

Author:

eugene

Message:

updated mkobj so it generates a resonable fake object

File:

• trunk/Ohana/src/opihi/mana/deimos_mkobj.c (modified) (9 diffs)

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

@@ -6 +6 @@
 
   // input parameters:
-  //   trace       : spline fit of slit central x pos vs y-coord
-  //   profile     : slit window profile (vector)
-  //   object      : vector of object flux vs y-coord
-  //   sky         : vector of local sky signal vs y-coord
-  //   background  : vector of extra-slit background flux vs y-coord
-  //   PSF         : point-spread function vector (flux normalized, x-dir)
-  //   LSF         : sline-spread function vector (flux normalized, y-dir)
-  //   stilt       : slit tilt response : 2D kernel? 
+  // * trace       : spline fit of slit central x pos vs y-coord
+  // * profile     : slit window profile (vector)
+  // * object      : vector of object flux vs y-coord 
+  // * sky         : vector of local sky signal vs y-coord
+  // * background  : vector of extra-slit background flux vs y-coord
+  // * PSF         : point-spread function vector (flux normalized, x-dir)
+  // * LSF         : sline-spread function vector (flux normalized, y-dir)
+  // * stilt       : slit tilt response : 2D kernel? 
 
   int N;
@@ -30 +30 @@
   Buffer *output  = NULL;
 
+  float stilt = 0.0; // angle of the slit
+
   if ((N = get_argument (argc, argv, "-trace"))) {
     remove_argument (N, &argc, argv);
@@ -63 +65 @@
     remove_argument (N, &argc, argv);
     if ((lsf = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-stilt"))) {
+    remove_argument (N, &argc, argv);
+    stilt = atof (argv[N]);
     remove_argument (N, &argc, argv);
   }
@@ -74 +81 @@
   }
 
-  if (argc != 4) {
+  if (argc != 3) {
     gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-Nwave N] [-object vector] [-sky vector] [-backgnd vector] [-trace spline] [-profile vector]\n");
     return FALSE;
@@ -110 +117 @@
 
   if ((output = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE);
-  int Nx = atoi(argv[2]);
+  int NxBase = atoi(argv[2]);
+  if (NxBase % 2 == 0) NxBase ++; // force Nx to be odd
+  int Nx = NxBase;
   
+  // if we are appying a trace offset spline, we need to generate an output window which
+  // is Nx + the full swing of the trace, then window back down
+  if (trace) {
+    float dXmin = +1000;
+    float dXmax = -1000;
+    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);
+      dXmin = MIN (dx, dXmin);
+      dXmax = MAX (dx, dXmax);
+    }
+    int dXrange = ceil(dXmax - dXmin);
+    if (dXrange % 2) dXrange ++; // force dXrange to be even (so Nx will be odd)
+    Nx += dXrange;
+  }
+
   ResetBuffer (output, Nx, Ny, -32, 0.0, 1.0);
 
@@ -117 +142 @@
   opihi_flt *skyV = (opihi_flt *) (sky    ? sky->elements.Flt : NULL);
   opihi_flt *psfV = (opihi_flt *) (psf    ? psf->elements.Flt : NULL);
-
-  float *out = (float *) output[0].matrix.buffer;
+  opihi_flt *lsfV = (opihi_flt *) (lsf    ? lsf->elements.Flt : NULL);
 
   // psf vector must have odd number of pixels:
   int Npsf = 0;
-  int Npof = floor(Nx/2);
+  int NpsfFull = 0;
   if (psf) {
     if (psf->Nelements % 2 == 0) {
@@ -128 +152 @@
       return FALSE;
     }
-    Npsf = psf->Nelements / 2;
-    Npof -= Npsf;
-  }
+    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);
 
   // loop over the y positions
@@ -138 +178 @@
     opihi_flt skyVy = skyV ? skyV[iy] : 0.0; // if sky is not supplied, flux defaults to 0.0
       
+    int Npof = floor(Nx/2);
+    if (psf) Npof -= Npsf;
+
     // flux = obj * PSF + sky
     for (int ix = 0; ix < Nx; ix++) {
@@ -144 +187 @@
 
       if (psfV) {
-        int n = ix - Npof; // n is the pixel in the PSF corresponding to the ix pixel in the output image
-        if (n < 0) continue;
-        if (n >= Npsf) continue;
-        
-        value += objVy * psfV[n];
+        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;
       }   
-      out[ix + iy*Nx] = value;
-    }
-  }
+      s1[ix + iy*Nx] = value;
+    }
+  }
+
+  // *** 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]);
+    }
+  }
+
+  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;
+  }
+
   return TRUE;
 }