baseline.c

Go to the documentation of this file.

/* 
Copyright 2021 Arizona Board of Regents on behalf of the University of Arizona.

For commercial uses, to obtain a license to sell and or sublicense copies of
the Software please contact the University of Arizona at Tech Launch Arizona:
lewish@tla.arizona.edu.

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction for non-commercial use educational and
research purposes, including without limitation the rights to use, copy,
modify, merge, publish, distribute, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:

The entire above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

1. Additional Required Provisions

1.1.  Arbitration. The parties agree that if a dispute arises between them
concerning this Agreement, the parties may be required to submit the matter to
arbitration pursuant to Arizona law.

1.2.  Applicable Law and Venue. This Agreement shall be interpreted pursuant
to the laws of the State of Arizona. Any arbitration or litigation between the
Parties shall be conducted in Pima County, ARIZONA, and LICENSEE hereby
submits to venue and jurisdiction in Pima County, ARIZONA.

1.3.  Non-Discrimination. The Parties agree to be bound by state and federal
laws and regulations governing equal opportunity and non-discrimination and
immigration.

1.4.   Appropriation of Funds. The Parties recognize that performance by
ARIZONA may depend upon appropriation of funds by the State Legislature of
ARIZONA. If the Legislature fails to appropriate the necessary funds, or if
ARIZONA’S appropriation is reduced during the fiscal year, ARIZONA may cancel
this Agreement without further duty or obligation. ARIZONA will notify
LICENSEE as soon as reasonably possible after it knows of the loss of funds.

1.5.   Conflict of Interest. This Agreement is subject to the provisions of
A.R.S. 38-511 and other conflict of interest regulations. Within three years
of the EFFECTIVE DATE, ARIZONA may cancel this Agreement if any person
significantly involved in initiating, negotiating, drafting, securing, or
creating this Agreement for or on behalf of ARIZONA becomes an employee or
consultant in any capacity of LICENSEE with respect to the subject matter of
this Agreement.

 */


#include "extern.h"

#define NUM_OF_BCHANNEL_FOR_POLY 10

int    valid_bl_averages;
double bl_averages[MAX_BASELINES];


int zero_bl_averages()
{
  bzero((char *)&bl_averages[0], sizeof(bl_averages));

  valid_bl_averages = 0;

  return(0);
}


int doBaseline()
{
  if(env->mode == MODE_CONT)
  {
    mapBaseline();
  }
  else
  {
    env->dobaseline = 1;
    env->bshow      = 0;

    if(isSet(PF_G_FLAG))
    {
      doGraceFit();
    }
    else
    {
      line_doBaseLine(h0);
    }
  }

  return(0);
}


int computeEachBaseLine(image)
struct LP_IMAGE_INFO *image;
{
  int    a, l, j, k, n, foundone=0;
  float  *p;

  zero_bl_averages();

  p = image->yplot1;

  n = getNumericFlag(); /* Does the user want to use just one of the regions for average? */

  if(n < 0) /* No, User hasn't asked for only one region; average them all */
  {
    a = 0;
    env->rst = 0.0;

    for(l=0;l<MAX_BASELINES;l++)
    {
      if(env->baselines[l][0] <= 0 || env->baselines[l][1] <= 0)
      {
        continue; /* Skip it */
      }

      if(env->baselines[l][0] > image->head.noint || env->baselines[l][1] > image->head.noint)
      {
        continue;
      }

      k = 0;

      for(j=env->baselines[l][0];j<env->baselines[l][1];j++)
      {
        if(j > env->bdrop && j < image->head.noint - env->edrop)
        {
          bl_averages[l] += p[j];
          k++;

          foundone++;
        }
      }

      if(k)
      {
        bl_averages[l] /= (double)k;

        env->rst += bl_averages[l];
        a++;
      }
    }

    if(a)
    {
      env->rst /= (double)a;
    }

    valid_bl_averages = 1;
  }
  else  /* Yes */
  {
    k = 0;

    for(j=env->baselines[n][0];j<env->baselines[n][1];j++)
    {
      if(j > env->bdrop && j < image->head.noint - env->edrop)
      {
        bl_averages[n] += p[j];
        k++;

        foundone++;
      }
    }

    if(k)
    {
      bl_averages[n] /= (double)k;

      env->rst = bl_averages[n];
    }
  }

  k = 0;
  for(l=0;l<MAX_BASELINES;l++)
  {
    if(bl_averages[l] != 0.0)
    {
      st_report("computeEachBaseLine(): Average values for Baseline Region: %d, %f", l, bl_averages[l]);
    }
  }

  return(foundone);
}



int computeBaseLine(image)
struct LP_IMAGE_INFO *image;
{
  int i, j, l, k, m=0;
  float      x[MAX_CHANS],
           sig[MAX_CHANS],
      lsqArray[MAX_CHANS];
  float a, b, siga1, sigb1, chi2_1, q1;
  float *p, r;
  double avg = 0.0;

  bzero((char *)&x,        sizeof(x));
  bzero((char *)&sig,      sizeof(sig));
  bzero((char *)&lsqArray, sizeof(lsqArray));

  p = image->yplot1;

  k=0;

  for(l=0;l<MAX_BASELINES;l++)
  {
    if(env->baselines[l][0] <= 0 || env->baselines[l][1] <= 0)
    {
      continue; /* Skip it */
    }

    if(env->baselines[l][0] > image->head.noint || env->baselines[l][1] > image->head.noint)
    {
      continue;
    }

    m++;
    for(j=env->baselines[l][0];j<env->baselines[l][1];j++)
    {
      if(j > env->bdrop && j < image->head.noint - env->edrop)
      {
        sig[k]      = 0.0;
        x[k]        = (float)j;
        lsqArray[k] = p[j];
        k++;
      }
    }
  }

  if(m == 0)
  {
    st_fail("No baseline regions defined");
    return(0);
  }

  for(i=0;i<k;i++)
  {
    avg += lsqArray[i];
  }

  avg /= (double)k;

  if(isSet(PF_N_FLAG))
  {
    env->rst = avg;

    return(k);
  }

  if(k>=4)
  {
    sptip_lsqfit(x, lsqArray, 0, k, sig, 0, &b, &a, &sigb1, &siga1, &chi2_1, &q1);

    computeRms(lsqArray, k, &r);

    image->binfo.a   = a;
    image->binfo.b   = b;
    image->binfo.rms = r;
  }
  else
  {
    image->binfo.a   = 0.0;
    image->binfo.b   = 0.0;
    image->binfo.rms = 100.0;
  }

  env->rms       = image->binfo.rms;
  env->slope     = image->binfo.a;
  env->intercept = image->binfo.b;

  return(k);
}



/* Take the yplot1 array in image and remove a baseline */
int line_doBaseLine(image)
struct LP_IMAGE_INFO *image;
{
  int i, ret, first = 1;
  float *f;

  copyEnv2Image(image); /* Make sure the structs are up to date */

  /* Compute baseline */
  ret = computeBaseLine(image);

  /* Remove any baseline */
  if(image->binfo.doBaseLine)
  {
    if(first && ret)
    {
      st_report("Removing Baselines from %s", image->name);
      first = 0;
    }

    if(verboseLevel & VERBOSE_BASELINES)
    {
      st_report("Slope %f, Intercept %f", image->binfo.a, image->binfo.b);
    }

    f = &image->yplot1[0];
    for(i=0;i<image->head.noint;i++,f++)
    {
      /* Do not remove if it's already zero; otherwise 
         you get bogus results from shifted data which
         zeros out the unshifted portion */
      if(*f != 0.0)
      {
        *f = *f - ((float)i * image->binfo.a + image->binfo.b);
      }
    }
  }

  return(0);
}


int isInBaseLineRegion(chan)
int chan;
{
  int i;

  for(i=0;i<MAX_BASELINES;i++)
  {
    if(env->baselines[i][0] <= 0 || env->baselines[i][1] <= 0)
    {
      continue; /* Skip it */
    }

    if(chan >= env->baselines[i][0] && chan <= env->baselines[i][1])
    {
      return(1);
    }
  }

  return(0);
}


int tellBaselines()
{
  int i;

//  st_report("Remove Baselines:               %s", env->dobaseline == 1 ? "Yes" : " No");

  for(i=0;i<MAX_BASELINES;i++)
  {
    if(i==0)
    {
      st_report("Baseline Regions:   [%2d]   %5d %5d", i, env->baselines[i][0], env->baselines[i][1]);
    }
    else
    {
      if(env->baselines[i][0]+env->baselines[i][1] > 0)
      {
        st_report("                    [%2d]   %5d %5d", i, env->baselines[i][0], env->baselines[i][1]);
      }
    }
  } 

  return(0);
}


int setNegBaselines()
{
  int i, n, c1 = 0, c2 = 0;
  int base, noint, old_bmark;
  int j=0, k=0;
  int  vals[2][2];
  double  arg1, arg2;
  char hint[256];
  struct LP_IMAGE_INFO *image = h0;

  noint      = image->head.noint;

  if(verboseLevel & VERBOSE_BASELINES)
  {
    st_report("Got %d tokens", args->ntokens);
  }

  if(args->ntokens > 2) /* Did user pass args */
  {
    arg1 = atof(args->secondArg);
    arg2 = atof(args->thirdArg);

    if(verboseLevel & VERBOSE_BASELINES)
    {
      st_report("Setting NBase to %.3f %.3f", arg1, arg2);
    }

    if(image->x1data[0] < image->x1data[noint-1])       /* X-Axis Increases */
    {
      if(arg1 > arg2)                                   /* Args in reverse order */
      {
        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("X increasing, args decreaing");
        }

        parseMouseValues(arg2, 0.0, 0.0, 0.0);

        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("Channel: %.0f", env->ccur);
        }

        c1 = (int)env->ccur;

        parseMouseValues(arg1, 0.0, 0.0, 0.0);

        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("Channel: %.0f", env->ccur);
        }

        c2 = (int)env->ccur;
      }
      else                                              /* Args in correct order */
      {
        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("X increasing, args increasing");
        }

        parseMouseValues(arg1, 0.0, 0.0, 0.0);

        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("Channel: %.0f", env->ccur);
        }

        c1 = (int)env->ccur;

        parseMouseValues(arg2, 0.0, 0.0, 0.0);

        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("Channel: %.0f", env->ccur);
        }

        c2 = (int)env->ccur;
      }

      if(c1 > -1 && c2 < image->head.noint-2)
      {
        bzero((char *)&env->baselines[0], sizeof(env->baselines));

        env->baselines[0][0] = 1;
        env->baselines[0][1] = c1;

        env->baselines[1][0] = c2;
        env->baselines[1][1] = image->head.noint-1;

        tellBaselines();

        return(0);
      }
      else
      {
        st_fail("Error in Translating Channels: %d %d", c1, c2);

        return(1);
      }
    }
    else                                                /* X-Axis decreases */
    {
      if(arg1 > arg2)                                   /* Args in correct order */
      {
        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("X decreasing, args decreaing");
        }

        parseMouseValues(arg1, 0.0, 0.0, 0.0);

        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("Channel: %.0f", env->ccur);
        }

        c1 = (int)env->ccur;

        parseMouseValues(arg2, 0.0, 0.0, 0.0);

        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("Channel: %.0f", env->ccur);
        }

        c2 = (int)env->ccur;
      }
      else                                              /* Args in reverse order */
      {
        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("X decreasing, args increasing");
        }

        parseMouseValues(arg2, 0.0, 0.0, 0.0);

        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("Channel: %.0f", env->ccur);
        }
        c1 = (int)env->ccur;

        parseMouseValues(arg1, 0.0, 0.0, 0.0);

        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("Channel: %.0f", env->ccur);
        }

        c2 = (int)env->ccur;
      }

      if(c1 > -1 && c2 < image->head.noint-2)
      {
        bzero((char *)&env->baselines[0], sizeof(env->baselines));

        env->baselines[0][0] = 1;
        env->baselines[0][1] = c1;

        env->baselines[1][0] = c2;
        env->baselines[1][1] = image->head.noint-1;

        tellBaselines();

        return(0);
      }
      else
      {
        st_fail("Error in Translating Channels: %d %d", c1, c2);

        return(1);
      }
    }
  }

  old_bmark  = env->bmark;
  env->bmark = 0;

  xx();

  n = 2;

  bzero((char *)&vals[0], sizeof(vals));

  for(i=0;i<n;i++)
  {
    sprintf(hint, "Click Left Mouse Button on %s Edge of Segment %d", !(i%2) ? "Left" : "Right", j+1);
    _doCCur(hint);

    base = (int)round(env->ccur);

    if(base < 1)
    {
      base = 1;
    }

    if(base >= noint)
    {
      base = noint-1;
    }

    vals[j][k] = base;
    k++;

    if(k == 2)
    {
      j++;
      k = 0;
    }
  }

  bzero((char *)&env->baselines[0], sizeof(env->baselines));

  env->baselines[0][0] = 1;
  env->baselines[0][1] = vals[0][0];

  env->baselines[1][0] = vals[0][1];
  env->baselines[1][1] = image->head.noint-1;

  tellBaselines();

  env->bmark = old_bmark;

  return(0);
}


int sortBaselines()
{
  int i, k, j=0, tmp;
  int bl[MAX_BASELINES][2];
  int avgs[MAX_BASELINES];
  struct LP_IMAGE_INFO *image = h0;

  bcopy((char *)&env->baselines[0][0], (char *)&bl[0][0], sizeof(bl));
  bzero((char *)&env->baselines[0][0], sizeof(env->baselines));
  bzero((char *)&avgs[0], sizeof(avgs));

  /* Find the average of each baseline segment */
  for(i=0;i<MAX_BASELINES;i++)
  {
    /* should be in numeric order Left to Right */
    if(bl[i][0] > bl[i][1])
    {
      st_report("Flipping baseline order for Segment %d", i+1);

      tmp = bl[i][0];

      bl[i][0] = bl[i][1];
      bl[i][1] = tmp;
    }

    avgs[i] = (bl[i][0] + bl[i][1]) / 2;
  }

  /* Now, for all possible channel numbers, place the segments in order */
  for(k=0;k<MAX_CHANS;k++)
  {
    if(k > image->head.noint) /* Can't be higher then this */
    {
      break;
    }

    for(i=0;i<MAX_BASELINES;i++)
    {
      if(avgs[i] == 0) /* Not set */
      {
        continue;
      }

      if(avgs[i] == k)
      {
        env->baselines[j][0] = bl[i][0];
        env->baselines[j][1] = bl[i][1];
        j++;
      }
    }
  }

  return(0);
}


int findEmptyBaseline()
{
  int i;

  for(i=0;i<MAX_BASELINES;i++)
  {
    if(env->baselines[i][0] + env->baselines[i][1] == 0)
    {
      return(i);
    }
  }

  return(-1);
}


int addOneOrMoreSegment()
{
  int i, j, n, new, ret, base, k=0;
  char hint[256];
  struct LP_IMAGE_INFO *image = h0;

  if(args->ntokens > 1) /* Passed param? */
  {
    n = atoi(args->secondArg);
  }
  else
  {
    n = 1;
  }

  for(j=0;j<n;j++)
  {
    new = findEmptyBaseline();

    if(new < 0)
    {
      st_fail("Baselines are full");

      sortBaselines();

      return(1);
    }

    for(i=0;i<2;i++)
    {
      sprintf(hint, "Click Left Mouse Button on %s Edge of Segment %d", !(i%2) ? "Left" : "Right", new+1);
      ret = _doCCur(hint);

      if(ret)
      {
        return(1);
      }

      base = (int)round(env->ccur);

      if(base < 1)
      {
        base = 1;
      }

      if(base >= image->head.noint)
      {
        base = image->head.noint-1;
      }

      if(base > -1)
      {
        env->baselines[new][k] = base;
        k++;
      }
    }
  }

  sortBaselines();

  return(0);
}



int removeOneSegment()
{
  int i, ret, now;
  char hint[256];

  sprintf(hint, "Click Left Mouse Button Somewhere in Segment to Remove");
  ret = _doCCur(hint);

  if(ret)
  {
    return(1);
  }

  now = (int)round(env->ccur);

  for(i=0;i<MAX_BASELINES;i++)
  {
    if(env->baselines[i][0] && env->baselines[i][1])
    {
      if(now >= env->baselines[i][0] && now <= env->baselines[i][1])
      {
        env->baselines[i][0] = 0;
        env->baselines[i][1] = 0;

        sortBaselines();

        return(0);
      }
    }
  }

  return(0);
}


int setBaselines()
{
  int i, n, base, j=0, k=0, ret, old_bmark;
  char hint[256];
  struct LP_IMAGE_INFO *image = h0;

  if(isSet(PF_C_FLAG)) /* Clear all and return */
  {
    st_report("Clearing out baselines segments");

    /* Zero out the current baselines */
    bzero((char *)&env->baselines[0], sizeof(env->baselines));

    return(0);
  }

  if(isSet(PF_A_FLAG)) /* Add a baseline region */
  {
    addOneOrMoreSegment();

    tellBaselines();

    return(0);
  }

  if(isSet(PF_R_FLAG)) /* Remove a baseline region */
  {
    removeOneSegment();

    tellBaselines();

    return(0);
  }

  if(args->ntokens > 1) /* Passed param? */
  {
    n = atoi(args->secondArg);
  }
  else
  {
    sprintf(hint, "Enter Number of Baseline Regions [ Up to %d]: ", MAX_BASELINES);

    do
    {
      n = getNumericResponse(hint);
    }
    while(n < 0 && n > MAX_BASELINES);
  }

  if(n == 0)
  {
    return(0);
  }

  if(n > MAX_BASELINES)
  {
    st_fail("Too Many Baseline Regions: [0 - %d]", MAX_BASELINES);

    return(1);
  }

  n *= 2;

  /* Zero out the current baselines */
  bzero((char *)&env->baselines[0], sizeof(env->baselines));

  if(n == 0)
  {
    st_print("Clearing out defined baseline regions\n");
    env->dobaseline = 0;

    return(0);
  }

  old_bmark = env->bmark;
  env->bmark = 0;
  xx();

  for(i=0;i<n;i++)
  {
    sprintf(hint, "Click Left Mouse Button on %s Edge of Segment %d", !(i%2) ? "Left" : "Right", j+1);
    ret = _doCCur(hint);

    if(ret)
    {
      sortBaselines();

      return(1);
    }

    base = (int)round(env->ccur);

    if(base < 1)
    {
      base = 1;
    }

    if(base >= image->head.noint)
    {
      base = image->head.noint-1;
    }

    if(base > -1)
    {
      env->baselines[j][k] = base;
      k++;

      if(k == 2)
      {
        j++;
        k = 0;
      }

      if(j == MAX_BASELINES)
      {
        if(old_bmark)
        {
          env->bmark = old_bmark;
        }
        else
        {
          env->bmark = 1;
        }

        sortBaselines();

        tellBaselines();

        return(0);
      }
    }
    else
    {
      if(old_bmark)
      {
        env->bmark = old_bmark;
      }
      else
      {
        env->bmark = 1;
      }

      st_fail("Can't set baseline channel of < 0");

      return(1);
    }
  }

  if(old_bmark)
  {
    env->bmark = old_bmark;
  }
  else
  {
    env->bmark = 1;
  }

  sortBaselines();

  tellBaselines();

  return(0);
}


int baselinesDefined(image)
struct LP_IMAGE_INFO *image;
{
  int i;

  for(i=0;i<MAX_BASELINES;i++)
  {
    if(image->binfo.baselines[i][0] == 0 && image->binfo.baselines[i][1] == 0)
    {
      continue;
    }
    else
    {
      /* found at least one segment defined */
      return(1);
    }
  }

  return(0);
}


int findLastGood(image, indx, lgn, lg)
struct LP_IMAGE_INFO *image;
int indx, *lgn;
float *lg;
{
  int i, n;
  float *f;

  n = 0;

  f = &image->bplot[indx];
  for(i=indx;i>=n;i--,f--)
  {
    if(*f != 0.0)
    {
      *lg = *f;
      *lgn = i;

      return(0);
    }
  }

  return(1);
}


int findNextGood(image, indx, ngn, ng)
struct LP_IMAGE_INFO *image;
int indx, *ngn;
float *ng;
{
  int i, n;
  float *f;


  n = image->head.noint;

  f = &image->bplot[indx];
  for(i=indx;i<n;i++,f++)
  {
    if(*f != 0.0)
    {
      *ng = *f;
      *ngn = i;

      return(0);
    }
  }

  return(0);
}


int doBshape()
{
  int b, n, b1, b2, n1, n2, i, j, k;
  int lgIndx, nxIndx, idiff, noint;
  int lastbl = 0, recording=0, start, end;
  char nstr[80], bname[256];
  float lg, ng, diff, unit;
  float lastGoodVal, firstGoodVal;
  double xd[MAX_CHANS], yd[MAX_CHANS], coeff[32];
  struct LP_IMAGE_INFO *image = h0;             /* h0 is a pointer to the current Hx */
  struct LP_IMAGE_INFO *b8image = &H8;          /* H8 the real thing */
  struct LP_IMAGE_INFO *b9image = &H9;          /* H9 the real thing */
  FILE *fp;

  bzero((char *)xd,    sizeof(xd));
  bzero((char *)yd,    sizeof(yd));
  bzero((char *)coeff, sizeof(coeff));
  bzero((char *)image->bplot, sizeof(image->bplot));

  if(!baselinesDefined(image))
  {
    st_fail("No baseline regions defined");
    return(1);
  }

  if(env->nbshape < 2)
  {
    st_fail("Error: nbshape must be greater than 1");
    return(1);
  }

  if(isSet(PF_R_FLAG))
  {
    recording = 1;
  }

  noint = (int)image->head.noint;

  st_report("Computing baseline using nbshape = %d", env->nbshape);

/* First perform a boxcar smoothing on the data withing each baselines regions */
  for(i=0;i<MAX_BASELINES;i++)
  {
    b = image->binfo.baselines[i][0];
    n = image->binfo.baselines[i][1] - image->binfo.baselines[i][0];

    if(n)
    {
      boxcar(&image->yplot1[b], &image->bplot[b], n, env->nbshape);

      lastbl = i;
    }
  }

  image->pinfo.bplot = 1;

    /* Now bplot array contains a baseline profile with areas not
       defined as baseline areas zeroed out. Think areas with lines.
       Now fill in the areas between baseline areas */

  if(env->bsmoothing == LINEAR_BASELINE_SMOOTHING)
  {
    st_report("Performing Linear Smoothing between Baseline Regions");

    for(i=0;i<noint;i++)
    {
      if(i == 0)
      {
        if(image->bplot[i] == 0.0)
        {
          findNextGood(image, i, &nxIndx, &ng);
          image->bplot[i] = ng;
        }
      }
      else
      if(i == noint-1)
      {
        if(image->bplot[i] == 0.0)
        {
          findLastGood(image, i, &lgIndx, &lg);
          image->bplot[i] = lg;
        }
      }
      else
      if(image->bplot[i] == 0.0)
      {
        findLastGood(image, i, &lgIndx, &lg);
        findNextGood(image, i, &nxIndx, &ng);

        diff  = lg - ng;
        idiff = lgIndx - nxIndx;

        if(idiff)
        {
          unit = diff / (float)idiff;

          if(verboseLevel & VERBOSE_BASELINES)
          {
            st_report("Missing Chan = %4d", i);
            st_report("Last good    = %4d, %f", lgIndx, lg);
            st_report("Next good    = %4d, %f", nxIndx, ng);
            st_report("Total Diff   = %f", diff);
            st_report("Unit Diff    = %f", unit);
          }

          image->bplot[i] = lg + unit * (double)(i-lgIndx);
        }
      }
    }
  }
  else
  if(env->bsmoothing == POLY_BASELINE_SMOOTHING)
  {
    if(!(env->nfit > 0 && env->nfit < MAX_NFIT))
    {
      st_fail("nfit must be > 0 & < %d; not %d", MAX_NFIT, env->nfit);
      return(1);
    }

    switch(env->nfit)
    {
      case  1: strcpy(nstr, "st"); break;
      case  2: strcpy(nstr, "nd"); break;
      case  3: strcpy(nstr, "rd"); break;

      default: 
        strcpy(nstr, "th"); 
        break;
    }

    st_report("Performing %d%s Order Polynomial Smoothing between Baseline Regions", env->nfit, nstr);

    /* Now fill in each region between baselines with a polynomial interp */
    for(i=0;i<(MAX_BASELINES);i++)
    {
      b1 = image->binfo.baselines[i][0];
      n1 = image->binfo.baselines[i][1] - image->binfo.baselines[i][0];

      b2 = image->binfo.baselines[i+1][0];
      n2 = image->binfo.baselines[i+1][1] - image->binfo.baselines[i+1][0];

      if(n1>4 && n2>4)
      {
        if(recording)
        {
          sprintf(bname, "bline-poly-%02d.dat", i);

          if((fp = fopen(bname, "w")) == NULL)
          {
            st_fail("Unable to open %s", bname);
            recording = 0;
          }
        }

        bzero((char *)xd,    sizeof(xd));
        bzero((char *)yd,    sizeof(yd));
        bzero((char *)coeff, sizeof(coeff));

        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("Interp between region %2d & %2d; %4d[%4d] to %4d[%4d]", i, i+1, b1, n1, b2, n2);
        }

        /* I only want NUM_OF_BCHANNEL_FOR_POLY channels in either side of the gap; if possible */

        j = 0;
        if(n1 > (NUM_OF_BCHANNEL_FOR_POLY-1) && n2 > (NUM_OF_BCHANNEL_FOR_POLY-1))      /* Both have plenty of channels */
        {
          start = b1 + n1 - NUM_OF_BCHANNEL_FOR_POLY;
          end   = b1 + n1;

          for(k=start;k<end;k++)
          {
            xd[j] = image->x1data[k];
            yd[j] = image->bplot[k];
            j++;
          }

          start = b2;
          end   = b2 + NUM_OF_BCHANNEL_FOR_POLY;

          for(k=start;k<end;k++)
          {
            xd[j] = image->x1data[k];
            yd[j] = image->bplot[k];
            j++;
          }
        }
        else
        if(n1 > (NUM_OF_BCHANNEL_FOR_POLY-1))                                           /* n2 somewhere between 2 - 9 channels */
        {
          start = b1 + n1 - NUM_OF_BCHANNEL_FOR_POLY;
          end   = b1 + n1;

          for(k=start;k<end;k++)
          {
            xd[j] = image->x1data[k];
            yd[j] = image->bplot[k];
            j++;
          }

          start = b2;
          end   = b2 + n2;

          for(k=start;k<end;k++)
          {
            xd[j] = image->x1data[k];
            yd[j] = image->bplot[k];
            j++;
          }
        }
        else
        if(n2 > (NUM_OF_BCHANNEL_FOR_POLY-1))                                           /* n1 somewhere between 2 - 9 channels */
        {
          start = b1;
          end   = b1 + n2;

          for(k=start;k<end;k++)
          {
            xd[j] = image->x1data[k];
            yd[j] = image->bplot[k];
            j++;
          }

          start = b2;
          end   = b2 + NUM_OF_BCHANNEL_FOR_POLY;

          for(k=start;k<end;k++)
          {
            xd[j] = image->x1data[k];
            yd[j] = image->bplot[k];
            j++;
          }
        }
        else                                                                    /* both n1 & n2 somewhere between 2 - 9 channels */
        {
          for(k=b1;k<(b1+n1);k++)
          {
            xd[j] = image->x1data[k];
            yd[j] = image->bplot[k];
            j++;
          }

          for(k=b2;k<(b2+n2);k++)
          {
            xd[j] = image->x1data[k];
            yd[j] = image->bplot[k];
            j++;
          }
        }

        if(j < 10)
        {
          st_fail("Not enought Baseline channels set for this region");
          continue;
        }

        fitcurve(xd, yd, j, env->nfit, coeff);

        /* Now fill in bplot with results */

        if(verboseLevel & VERBOSE_BASELINES)
        {
          st_report("%02d: Filling in channels %d to %d", i, b1+n1, b2);
        }

        for(k=(b1+n1-1);k<(b2+1);k++)
        {
          image->bplot[k] = leasev(coeff, env->nfit, image->x1data[k]);
        }

        if(recording)
        {
          for(k=b1;k<(b2+n2);k++)
          {
            fprintf(fp, "%f %f\n", image->x1data[k], leasev(coeff, env->nfit, image->x1data[k]));
          }

          fclose(fp);
        }

        /* Now if this is the first region, does it cover the first channels */
        if(i == 0 && b1 > 0)
        {
          if(env->bleadtrail)  /* Fill in leading channels with Polynomial from first good baseline ??? */
          {
            st_report("Poly-Filling in leading channels %d to %d", 0, b1);

            /* Recomute using entire 1st baseline region */
            j = 0;
            for(k=b1;k<b1+n1;k++)
            {
              xd[j] = image->x1data[k];
              yd[j] = image->bplot[k];
              j++;
            }

            fitcurve(xd, yd, j, env->nfit, coeff);

            j = 0;
            for(k=0;k<b1;k++)
            {
              image->bplot[k] = leasev(coeff, env->nfit, image->x1data[k]);
              j++;
            }
          }
          else /* Just take first good value and propagate to starting channel */
          {
            st_report("Linear-Filling in leading channels %d to %d", 0, b1);

            firstGoodVal = image->bplot[b1];

            j = 0;
            for(k=0;k<b1;k++)
            {
              image->bplot[k] = firstGoodVal;
              j++;
            }
          }

          if(verboseLevel & VERBOSE_BASELINES)
          {
            st_report("%02d: Filled in %d leading channels", i, j);
          }
        }

        /* Now, fill in trailing channels */

        if((i+1) == lastbl && (b2+n2) < noint)
        {
          lastGoodVal = image->bplot[b2+n2-1];

          if(env->bleadtrail) /* Continue Polynomial fitting out to the last channel ??? */
          {
            st_report("Poly-Filling in trailing channels %d to %d", b2+n2, noint);

          /* Recomute using entire 1st baseline region */
            j = 0;
            for(k=b2;k<b2+n2;k++)
            {
              xd[j] = image->x1data[k];
              yd[j] = image->bplot[k];
              j++;
            }

            fitcurve(xd, yd, j, env->nfit, coeff);

            j=0;
            for(k=b2+n2;k<noint;k++)
            {
              image->bplot[k] = leasev(coeff, env->nfit, image->x1data[k]);
              j++;
            }

            if(verboseLevel & VERBOSE_BASELINES)
            {
              st_report("Filled in %d trailing channels", j);
            }
          }
          else /* Take last good value and run it out to end */
          {
            st_report("Linear-Filling in trailing channels %d to %d", b2+n2, noint);

            j=0;
            for(k=b2+n2;k<noint;k++)
            {
              image->bplot[k] = lastGoodVal;
              j++;
            }
          }

          if(verboseLevel & VERBOSE_BASELINES)
          {
            st_report("Filled in %d trailing channels", j);
          }

        }
      }
      else
      {
        if(verboseLevel & VERBOSE_BASELINES)
        {
          if(n1 || n2)
          {
            st_report("Skipping short baseline region %d: n1 = %d, n2 = %d", i, n1, n2);
          }
        }
      }
    }
  }
  else
  {
    st_fail("bsmoothing can only be 0 or 1; not %d", env->bsmoothing);
    return(1);
  }

  if(isSet(PF_D_FLAG)) /* Duplicate to H8 image ?? */
  {
    st_report("Saving baseline model in H8");
    bcopy((char *)image, (char *)b8image, sizeof(struct LP_IMAGE_INFO));
    bcopy((char *)&image->bplot, (char *)b8image->yplot1, sizeof(b8image->yplot1));
    strcpy(b8image->name, "H8"); /* Opps, fix the name */
  }

  if(isSet(PF_S_FLAG))
  {
    args->pflags[PF_S_FLAG].set = 1;

    doBspline();

    if(isSet(PF_D_FLAG)) /* Duplicate to H9 image ?? */
    {
      st_report("Saving smoothed baseline model in H9");
      bcopy((char *)image, (char *)b9image, sizeof(struct LP_IMAGE_INFO));
      bcopy((char *)&image->bplot, (char *)b9image->yplot1, sizeof(b9image->yplot1));
      strcpy(b9image->name, "H9"); /* Opps, fix the name */
    }
  }

  return(0);
}


int do_contGraceFit()
{
  int    i, j, k, n, foundone = 0;
  float  *f;
  double xd[MAX_CHANS], yd[MAX_CHANS], coeff[32], avg;
  struct ACCUM_SCAN_RESULT  *p;
  struct ACCUM_RESULTS_SEND *s;

  if(isSet(PF_O_FLAG))          /* Perform this on the OBJECT arrays ?? */
  {
    return(do_objectGraceFit());
  }

  bzero((char *)xd,    sizeof(xd));
  bzero((char *)yd,    sizeof(yd));
  bzero((char *)coeff, sizeof(coeff));

  if(env->nfit == 0)
  {
    st_report("do_contGraceFit(): Using nfit = %d", env->nfit);

    for(j=0;j<MAX_IFS;j++)
    {
      n = accum_results_send.ar[j].num;

      if(n)
      {
        st_report("do_contGraceFit(): IF: %d contains %d records", j, n);

        p = &accum_results_send.ar[j].results[0];

        avg = 0.0;
        for(i=0;i<n;i++,p++)
        {
          switch(env->acc_plot_select_y)
          {
            case ACCUM_PLOT_Y_SIG:    avg += p->sig;        break;
            case ACCUM_PLOT_Y_SKY:    avg += p->sky;        break;
            case ACCUM_PLOT_Y_SMROR:  avg += p->smror;      break;
            case ACCUM_PLOT_Y_TSYS:   avg += p->tsys;       break;
            case ACCUM_PLOT_Y_TASTAR: avg += p->tastar;     break;
            case ACCUM_PLOT_Y_TR:     avg += p->tr;         break;
            case ACCUM_PLOT_Y_JANSKY: avg += p->janskys;    break;
            case ACCUM_PLOT_Y_NORMAL: avg += p->normalized; break;
          }
        }

        avg /= (double)n;

        st_report("do_contGraceFit(): Average for IF: (%d):, %f", j, avg);

        f = &accum_results_send.fitcurve[j][0];
        for(i=0;i<n;i++,f++)
        {
          *f = (float)avg;
        }

        foundone++;
      }
    }

    if(foundone)
    {
      accum_results_send.bplot  = 1;
      accum_results_send.ncoeff = env->nfit;

      sendAccumResults();
    }
    else
    {
      st_fail("do_contGraceFit(): No data to fit");
      return(1);
    }
  }
  else
  if(env->nfit > 0 && env->nfit < MAX_NFIT)
  {
    st_report("do_contGraceFit(): Using nfit = %d", env->nfit);

    /* Convert floats to doubles */

    for(j=0;j<MAX_IFS;j++)
    {
      n = accum_results_send.ar[j].num;

      if(n)
      {
        st_report("do_contGraceFit(): IF: %d contains %d records", j, n);

        p = &accum_results_send.ar[j].results[0];

        for(i=0;i<n;i++,p++)
        {
          switch(env->acc_plot_select_x)
          {
            case ACCUM_PLOT_X_EL:    xd[i] = p->el;     break;
            case ACCUM_PLOT_X_UT:    xd[i] = p->ut;     break;
            case ACCUM_PLOT_X_SCAN:  xd[i] = p->scan;   break;
            default:                 xd[i] = (double)i; break;
          }

          switch(env->acc_plot_select_y)
          {
            case ACCUM_PLOT_Y_SIG:    yd[i] = p->sig;        break;
            case ACCUM_PLOT_Y_SKY:    yd[i] = p->sky;        break;
            case ACCUM_PLOT_Y_SMROR:  yd[i] = p->smror;      break;
            case ACCUM_PLOT_Y_TSYS:   yd[i] = p->tsys;       break;
            case ACCUM_PLOT_Y_TASTAR: yd[i] = p->tastar;     break;
            case ACCUM_PLOT_Y_TR:     yd[i] = p->tr;         break;
            case ACCUM_PLOT_Y_JANSKY: yd[i] = p->janskys;    break;
            case ACCUM_PLOT_Y_NORMAL: yd[i] = p->normalized; break;
          }
        }

        fitcurve(xd, yd, n, env->nfit, coeff);

        if(env->issmt)
        {
          st_print("Linuxpops(C) >>  # do_contGraceFit(): Results for IF: (%s): ", mg_sideband[j]);
        }
        else
        {
          st_print("Linuxpops(C) >>  # do_contGraceFit(): Results for IF: (%s): ", kp_sideband[j]);
        }

        for(k=0;k<env->nfit+1;k++)
        {
          st_print("a%d = %.4e ", k, coeff[k]);

          accum_results_send.coeff[j][k] = coeff[k];
        }

        st_print("\n");

        /* Now fill in fitcurve with results */
        p = &accum_results_send.ar[j].results[0];
        f = &accum_results_send.fitcurve[j][0];

        for(i=0;i<n;i++,f++,p++)
        {
          *f = (float)leasev(coeff, env->nfit, p->el);
        }

        foundone++;
      }
    }

    if(foundone)
    {
      s = &accum_results_send;

      s->bplot = 1;
      s->ncoeff = env->nfit;

      sendAccumResults();

      if(isSet(PF_S_FLAG))              /* Save coeff to elgain terms?? */
      {
        st_report("Saving %d Elevation Gain Terms", env->nfit+1);

        env->elgain_term1 = 0.0;                /* Zero old ones out first */
        env->elgain_term2 = 0.0;
        env->elgain_term2 = 0.0;
        env->elgain_term4 = 0.0;
        env->elgain_term5 = 0.0;

        for(k=0;k<env->nfit+1;k++)
        {
          switch(k)
          {
            case 0: env->elgain_term1 = (s->coeff[0][k] + s->coeff[1][k] + s->coeff[2][k] + s->coeff[3][k]) / (float)foundone; break;
            case 1: env->elgain_term2 = (s->coeff[0][k] + s->coeff[1][k] + s->coeff[2][k] + s->coeff[3][k]) / (float)foundone; break;
            case 2: env->elgain_term3 = (s->coeff[0][k] + s->coeff[1][k] + s->coeff[2][k] + s->coeff[3][k]) / (float)foundone; break;
            case 3: env->elgain_term4 = (s->coeff[0][k] + s->coeff[1][k] + s->coeff[2][k] + s->coeff[3][k]) / (float)foundone; break;
            case 4: env->elgain_term5 = (s->coeff[0][k] + s->coeff[1][k] + s->coeff[2][k] + s->coeff[3][k]) / (float)foundone; break;
          }
        }
      }
    }
    else
    {
      st_fail("do_contGraceFit(): No data to fit");
      return(1);
    }
  }
  else
  {
    st_fail("Error: nfit out of range: %d, [0 - %d]", env->nfit, MAX_NFIT-1);
    return(1);
  }

  return(0);
}


int doGraceFit()
{
  int i, j, n, ret;
  struct LP_IMAGE_INFO *image = h0;
  double xd[MAX_CHANS], yd[MAX_CHANS], coeff[32], avg;

  if(env->mode == MODE_CONT)
  {
    ret = do_contGraceFit();

    return(ret);
  }

  bzero((char *)xd,    sizeof(xd));
  bzero((char *)yd,    sizeof(yd));
  bzero((char *)coeff, sizeof(coeff));

  n = image->head.noint;

  if(env->nfit == 0)
  {
    st_report("graceFit(): Using nfit = %d", env->nfit);

    image->pinfo.bplot = 1;

    j = 0;
    avg = 0.0;
    for(i=0;i<n;i++)
    {
      if(isInBaseLineRegion(i))
      {
        avg += image->yplot1[i];
        j++;
      }
    }

    if(j)
    {
      avg /= (double)j;

    /* Now fill in bplot with results */
      for(i=0;i<image->head.noint;i++)
      {
        image->bplot[i] = avg;
      }
    }
    else
    {
      st_fail("No baseline regions set");
      return(1);
    }
  }
  else
  if(env->nfit > 0 && env->nfit < MAX_NFIT)
  {
    st_report("graceFit(): Using nfit = %d", env->nfit);

    image->pinfo.bplot = 1;

    /* Convert floats to doubles */
    j = 0;
    for(i=0;i<n;i++)
    {
      if(isInBaseLineRegion(i))
      {
        xd[j] = image->x1data[i];
        yd[j] = image->yplot1[i];
        j++;
      }
    }

    if(j < 10)
    {
      st_fail("Not enought Baseline channels set");
      return(1);
    }

    fitcurve(xd, yd, j, env->nfit, coeff);

    /* Now fill in bplot with results */
    for(i=0;i<image->head.noint;i++)
    {
      image->bplot[i] = leasev(coeff, env->nfit, image->x1data[i]);
    }
  }
  else
  {
    st_fail("Error: nfit out of range: %d, [0 - 19]", env->nfit);
  }

  xx();

  return(0);
}


int doRmBase()
{
  double diff;
  struct LP_IMAGE_INFO *image = h0;

  st_report("Removing Computed Baseline from: %s", image->name);

  env->rm_fit_data = 1;

  removeFittedData(image);

  env->rm_fit_data = 0;

  if(isSet(PF_A_FLAG))
  {
    doAutoScale();
  }

  if(isSet(PF_F_FLAG))
  {
    diff = env->yupper - env->ylower;

    env->yupper =  (diff * 0.5);
    env->ylower = -(diff * 0.5);
  }

  return(0);
}


int doBmove()
{
  int i, ret, seg=0, now, then, min = 999999, mindx = -1;
  char hint[256];
  struct LP_IMAGE_INFO *image = h0;

  sprintf(hint, "Click Left Mouse Button on Edge of Segment to move");
  ret = _doCCur(hint);

  if(ret)
  {
    return(1);
  }

  now = (int)round(env->ccur);

  sprintf(hint, "Click Left Mouse Button in New Location");
  ret = _doCCur(hint);

  if(ret)
  {
    return(1);
  }

  then = (int)round(env->ccur);

  if(then < 1)
  {
    then = 1;
  }
  else
  if(then > image->head.noint)
  {
    then = image->head.noint-1;
  }

  for(i=0;i<MAX_BASELINES;i++)
  {
    if(env->baselines[i][0] != 0 && env->baselines[i][1] != 0) /* definded ?? */
    {
      if(abs(env->baselines[i][0]-now) < min)
      {
        seg = i;
        mindx = 0;
        min = abs(env->baselines[i][0]-now);
      }

      if(abs(env->baselines[i][1]-now) < min)
      {
        seg = i;
        mindx = 1;
        min = abs(env->baselines[i][1]-now);
      }
    }
  }

  if(seg >= 0 && seg < MAX_BASELINES)
  {
    if(mindx == 0 || mindx == 1)
    {
      st_report("Closest segment: %d, part: %d; Channel: %d", seg+1, mindx, env->baselines[seg][mindx]);
      st_report("Moving Segment %d edge from channel %d to %d", seg+1, now, then);

      env->baselines[seg][mindx] = then;

      if(isSet(PF_B_FLAG))
      {
        doBshape();
      }
      else
      if(isSet(PF_G_FLAG))
      {
        doGraceFit();
      }

      return(0);
    }
  }

  st_fail("Seg: %d or index: %d; out of range", seg, mindx);

  return(1);
}



int doBloop()
{
  while(1)
  {
    if(doBmove() == 1)
    {
      return(0);
    }

    if(isSet(PF_B_FLAG))
    {
      doBshape();
    }
    else
    if(isSet(PF_G_FLAG))
    {
      doGraceFit();
    }

    xx();
  }

  return(0);
}


int doNBaseset()
{
  int seg, start, end;

  /* /c flag clears the whole thing */
  if(isSet(PF_C_FLAG))
  {
    bzero((char *)&env->baselines[0], sizeof(env->baselines));
  }

  if(args->ntokens > 3)
  {
    seg = atoi(args->secondArg);

    if(seg >=0 && seg < MAX_BASELINES)
    {
      start = atoi(args->thirdArg);
      end   = atoi(args->fourthArg);

      st_report("Setting Seg %2d: %4d %4d", seg, start, end);

      env->baselines[seg][0] = start;
      env->baselines[seg][1] = end;
    }
    else
    {
      st_fail("Baseline Segment %d, out of range", seg);
      return(1);
    }
  }

  if(isSet(PF_L_FLAG))
  {
    tellBaselines();
  }

  return(0);
}


int writeBaseLines()
{
  int i, n=0;
  FILE *fp;

  if(args->ntokens < 2)
  {
    st_fail("Usage: writebase filename");
    return(0);
  }
  else
  {
    if((fp = fopen(args->secondArg, "w")) == NULL)
    {
      st_fail("Unable to open: %s", args->secondArg);
      return(0);
    }
    
    fprintf(fp, "nbaseset /c\n"); /* Clear them out first */

    for(i=0;i<MAX_BASELINES;i++)
    {
      if(env->baselines[i][0] ||  env->baselines[i][1])
      {
        fprintf(fp, "nbaseset %2d %4d %4d\n", i, env->baselines[i][0], env->baselines[i][1]);
        n++;
      }
    }

    st_report("Wrote %d baseline segments to: '%s'", n, args->secondArg);

    fclose(fp);
  }

  return(0);
}


int dcOffset()
{
  int    i, n, j;
  double avg;

  if(isSet(PF_A_FLAG)) /* take average and remove */
  {
    avg = 0.0;
    n   =  h0->head.noint;

    j = 0;
    for(i=0;i<n;i++)
    {
      if(isInBaseLineRegion(i))
      {
        avg += h0->yplot1[i];
        j++;
      }
    }

    if(j)
    {
      avg /= (double)j;
    }
    else
    {
      st_fail("No baseline regions defined");
      return(1);
    }

    for(i=0;i<n;i++)
    {
      h0->yplot1[i] -= avg;
    }

    return(0);
  }

  if(args->ntokens > 1)
  {
    avg = atof(args->secondArg);

    if(isnormal(avg))
    {
      n =  h0->head.noint;

      for(i=0;i<n;i++)
      {
        h0->yplot1[i] -= avg;
      }
    }
    else
    {
      printH_FlagHelp("dcoffset");

      return(1);
    }
  }
  else
  {
    printH_FlagHelp("dcoffset");

    return(1);
  }

  return(0);
}