IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 40755


Ignore:
Timestamp:
May 27, 2019, 8:42:23 AM (7 years ago)
Author:
eugene
Message:

coding trace analysis to handle bad edge cases; add red & blu sides for trace spline; adding variance to flat-field

Location:
branches/eam_branches/ohana.20190329/src/opihi
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r40662 r40755  
    2222  int N;
    2323  Buffer *image;
    24   Buffer *sig = NULL;
    25 
    26   if ((N = get_argument (argc, argv, "-sigma"))) {
    27     remove_argument (N, &argc, argv);
    28     if ((sig = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     24  Buffer *var = NULL;
     25
     26  if ((N = get_argument (argc, argv, "-variance"))) {
     27    remove_argument (N, &argc, argv);
     28    if ((var = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    2929    remove_argument (N, &argc, argv);
    3030  }
    3131
    3232  if (argc != 3) {
    33     gprint (GP_ERR, "USAGE: medimage add (name) (image) [-sigma sigma]\n");
     33    gprint (GP_ERR, "USAGE: medimage add (name) (image) [-variance variance]\n");
    3434    gprint (GP_ERR, "       add the given image to the set of images to be medianed\n");
    3535    gprint (GP_ERR, "       optionally supply variance image (for weighted calculations)\n");
    … …  
    3939  if ((image = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    4040
    41   if (sig) {
    42     if ((sig->matrix.Naxis[0] != image->matrix.Naxis[0]) ||
    43         (sig->matrix.Naxis[1] != image->matrix.Naxis[1])) {
    44       gprint (GP_ERR, "sigma buffer does not match image buffer dimensions\n");
     41  if (var) {
     42    if ((var->matrix.Naxis[0] != image->matrix.Naxis[0]) ||
     43        (var->matrix.Naxis[1] != image->matrix.Naxis[1])) {
     44      gprint (GP_ERR, "variance buffer does not match image buffer dimensions\n");
    4545      return FALSE;
    4646    }
    … …  
    6868  // new image should match existing medimage dimensions
    6969
    70   // AddMedImage (median, image, sig);
     70  // AddMedImage (median, image, var);
    7171  int Ninput = median->Ninput;
    7272  median->Ninput ++;
    7373  REALLOCATE (median->flx, float *, median->Ninput);
    74   REALLOCATE (median->sig, float *, median->Ninput);
     74  REALLOCATE (median->var, float *, median->Ninput);
    7575
    7676  ALLOCATE (median->flx[Ninput], float, median->Nx*median->Ny);
    7777  memcpy (median->flx[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
    7878
    79   median->sig[Ninput] = NULL;
    80   if (sig) {
    81     ALLOCATE (median->sig[Ninput], float, median->Nx*median->Ny);
    82     memcpy (median->sig[Ninput], sig->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
     79  median->var[Ninput] = NULL;
     80  if (var) {
     81    ALLOCATE (median->var[Ninput], float, median->Nx*median->Ny);
     82    memcpy (median->var[Ninput], var->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
    8383  }
    8484
    … …  
    105105    if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -irls, -wtmean\n"); return FALSE; }
    106106    mode = CALC_WTMEAN;
     107    remove_argument (N, &argc, argv);
     108  }
     109
     110  Buffer *variance = NULL;
     111  if ((N = get_argument (argc, argv, "-variance"))) {
     112    remove_argument (N, &argc, argv);
     113    if ((variance = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE);
    107114    remove_argument (N, &argc, argv);
    108115  }
    … …  
    145152        val[N] = v;
    146153        wgt[N] = 1.0;
    147         if (median->sig[n]) {
    148           float s = median->sig[n][Npix];
     154        if (median->var[n]) {
     155          float s = median->var[n][Npix];
    149156          if (!isfinite(s)) continue;
    150157          if (fabs(s) < 2*FLT_MIN) s = 2*FLT_MIN;
    151           wgt[N] = 1.0 / SQ(s);
     158          wgt[N] = 1.0 / s;
    152159        }
    153160        N++;

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

    r40662 r40755  
    4949  int Ny;
    5050  float **flx;
    51   float **sig;
     51  float **var;
    5252} MedImageType;
    5353

  • branches/eam_branches/ohana.20190329/src/opihi/lib.data/MedImageOps.c

    r40662 r40755  
    3333  for (i = 0; i < medimage[0].Ninput; i++) {
    3434    free (medimage[0].flx[i]);
    35     FREE (medimage[0].sig[i]);
     35    FREE (medimage[0].var[i]);
    3636  }
    3737  free (medimage[0].flx);
    38   free (medimage[0].sig);
     38  free (medimage[0].var);
    3939  free (medimage);
    4040}
    … …  
    7878  medimage->Ny = Ny;
    7979  ALLOCATE (medimage->flx, float *, 1);
    80   ALLOCATE (medimage->sig, float *, 1);
     80  ALLOCATE (medimage->var, float *, 1);
    8181
    8282  medimages[N] = medimage;

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

    r40754 r40755  
    7878 
    7979  float Sky_V = 0.0, BckL_V = 0.0, BckR_V = 0.0;
    80   int Sky_N = 0, BckL_N = 0, BckR_N = 0;
     80  float Sky_N = 0.0, BckL_N = 0.0, BckR_N = 0.0;
    8181  for (int i = 0; i < xprofile->Nelements; i++) {
     82    if (!isfinite(wprofile->elements.Flt[i])) continue;
     83    if (!isfinite(fprofile->elements.Flt[i])) continue;
     84    float wt = 1 / SQ(wprofile->elements.Flt[i]);
    8285    if (xprofile->elements.Flt[i] < Xs) {
    83       BckL_V += fprofile->elements.Flt[i];
    84       BckL_N ++;
     86      BckL_V += fprofile->elements.Flt[i]*wt;
     87      BckL_N += wt;
    8588      continue;
    8689    }
    8790    if (xprofile->elements.Flt[i] > Xe) {
    88       BckR_V += fprofile->elements.Flt[i];
    89       BckR_N ++;
     91      BckR_V += fprofile->elements.Flt[i]*wt;
     92      BckR_N += wt;
    9093      continue;
    9194    }
    92     Sky_V += fprofile->elements.Flt[i];
    93     Sky_N ++;
     95    Sky_V += fprofile->elements.Flt[i]*wt;
     96    Sky_N += wt;
    9497  }
    9598  float Sky = Sky_V / Sky_N;
    … …  
    211214  int Npts = 0;
    212215  for (int i = 0; i < NptsAll; i++) {
     216    if (!isfinite(wprofile[i])) continue;
     217    if (!isfinite(fprofile[i])) continue;
    213218    int keep = (mode == MODE_UPPER) ^ (xprofile[i] < limit);
    214219    if (keep) {

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

    r40659 r40755  
    77  // input parameters:
    88  //   slit profile response : vector of fractional flux vs x-coord
    9   //   trace                 : spline fit of slit central x pos vs y-coord
    10   //   flux                  : vector of input signal vs y-coord
     9  //   trace_ref,trace_blu   : spline fit of slit central x pos vs y-coord
     10  //   flux (optional)       : vector of input signal vs y-coord
     11  //   sky  (optional)       : vector of input signal vs y-coord
     12  //   fluxbins (optional)   : size of output image in y-direction
     13  //   NOTE: one of flux, sky, or fluxbins must be defined
     14
     15  // output : buff (an image with Nx defined by profile and Ny defined
     16  //                by flux, sky, or fluxbins)
    1117
    1218  int N;
    1319
    14   Vector *profile = NULL;
    15   Vector *flux    = NULL;
    16   Vector *sky     = NULL;
    17   Spline *trace   = NULL;
    18   Buffer *buff    = NULL;
     20  Vector *profile   = NULL;
     21  Vector *flux      = NULL;
     22  Vector *sky       = NULL;
     23  Spline *trace_red = NULL;
     24  Spline *trace_blu = NULL;
     25  Buffer *buff      = NULL;
    1926
    2027  // user-specified flux value (otherwise assumed to be 1)
    … …  
    4047  }
    4148
    42   if (argc != 4) {
    43     gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
     49  // for a red vs blu spline, we need to specify the split point
     50  // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
     51  // this more generic (need to define the ranges somewhere)
     52  int redlimit = 4096;
     53  if ((N = get_argument (argc, argv, "-redlimit"))) {
     54    remove_argument (N, &argc, argv);
     55    redlimit = atoi (argv[N]);
     56    remove_argument (N, &argc, argv);
     57  }
     58
     59  if (argc != 5) {
     60    gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace_red) (trace_blu) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
    4461    return FALSE;
    4562  }
    … …  
    5168
    5269  // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
    53   if ((profile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
    54   if ((trace   = FindSpline (argv[2])) == NULL) return (FALSE);
    55   if ((buff    = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE);
     70  if ((profile   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     71  if ((trace_red = FindSpline (argv[2])) == NULL) return (FALSE);
     72  if ((trace_blu = FindSpline (argv[3])) == NULL) return (FALSE);
     73  if ((buff      = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
    5674
    5775  // define the output window
    … …  
    7391
    7492    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
     93    Spline *trace = (iy < redlimit) ? trace_red : trace_blu;
    7594    float dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
    7695   
Note: See TracChangeset for help on using the changeset viewer.