 /*
  * Khoros: $Id: ldpeakpick.c,v 1.2 1992/03/20 23:28:54 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: ldpeakpick.c,v 1.2 1992/03/20 23:28:54 dkhoros Exp $";
#endif

 /*
  * $Log: ldpeakpick.c,v $
 * Revision 1.2  1992/03/20  23:28:54  dkhoros
 * VirtualPatch5
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as to the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including, for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: ldpeakpick.c
 >>>>
 >>>>      Program Name: dpeakpick
 >>>>
 >>>> Date Last Updated: Sat Feb 15 16:28:25 1992 
 >>>>
 >>>>          Routines: ldpeakpick - the library call for dpeakpick
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
struct p_info
{
  float p_value;
  int   p_position;
};

struct v_info
{
  float p2p_dist_mean,
        p2p_dist_var;
  int   p_count;
};

/* -library_includes_end */


/****************************************************************
*
* Routine Name: ldpeakpick - library call for dpeakpick
*
* Purpose:
*    
*    Find peaks or troughs in 1D data set
*    
*    

* Input:
*    
*    image(struct xvimage) - input data structure
*    
*    peak_flg
*         (int) - flag set to 1 if peaks are to be detected. Set to  0
*         if troughs are to be detected
*    
*    all_flg
*         (int) - flag set to 1 if all peaks are to be  detected.  Set
*         to 0 if user specifies the number of peaks to be detected.
*    
*    num_peaks
*         (int) - number of peaks/troughs to  detect  (only  valid  if
*         all_flg = 0)
*    
*    max_flg
*         (int) - if flag = 0 -> find first n peaks, if flag  =  1  ->
*         find max/min n peaks/troughs (only valid if all_flg = 0)
*    
*    end_flg
*         (int) - flag which signals whether the first and last points
*         in  a  data sequence should be considered for peak or trough
*         designation. If end_flg is 1, endpoints will  be  considered
*         for peak or trough designation, if 0, they will not.
*    
*    w_size(int) - number of steps on both sides of a point that  must
*         be monotonically increasing/decreasing before the point will
*         be counted as a peak or trough.
*    
*    hwr  (double) - minimum height to width ratio.  Point of interest
*         is  a  peak only if absolute value of the average slope over
*         window size is greater than hwr.
*    
*    process_dir
*         (int) - Direction of data processing: down  vectors  =  0  =
*         DSP_VECTOR, across bands = 1 = DSP_BAND
*    
*    

* Output:
*    
*    image(struct xvimage) - output peak/trough  implicit  VIFF  file.
*         All values are zero except peak/trough values.
*    
*    peak_info
*         (struct p_info) - 2D array containing peak/trough values and
*         positions.  If the data is complex, the number of vectors in
*         this array will be 2 times the  actual  number  of  vectors.
*         They  will  be  arranged  in  the  following order vect0:r,i
*         vect1:r,i ....
*    
*    vect_info
*         (struct v_info) - 1D array containing number of  peaks,  and
*         mean and variance of the distance between peaks/troughs, for
*         each vector.  If the data is complex, the number of  vectors
*         represented  in this array will be 2 times the actual number
*         of vectors.  The statistics will be arranged in the  follow-
*         ing order vect0:r,i vect1:r,i ....
*    
*    

*
* Written By: Donna Koechner
*    
*    

****************************************************************/


/* -library_def */
int
ldpeakpick(image, peak_flg, all_flg, num_peaks, max_flg, end_flg, w_size, 
           hwr, process_dir, peak_info, vect_info)
struct xvimage  *image;         /* input image structure */
int             peak_flg,       /* flg = 1 for peaks, 0 for troughs */
                all_flg,        /* flg = 1 find all peaks */
                num_peaks,      /* number of peaks or troughs to detect */
                max_flg,        /* flag = 0 -> find first n peaks,
                                   flag = 1 -> find max/min n peaks/troughs */
                end_flg,        /* if TRUE, end points can be peaks */
                w_size,         /* window size to each side of point */
                process_dir;    /* 0 = DSP_VECTOR, 1 = DSP_BAND */
double          hwr;            /* minimum height to width ratio */
struct p_info   ***peak_info;    /* 2D array of peak/trough values & positions */
struct v_info   **vect_info;     /* 1D array of vector peak statistics */


/* -library_def_end */

/* -library_code */
{
  char          **dload_vector(),
                *program = "ldpeakpick";
  int           num_vects,      /* number vectors returned from dloadvector */
                l_vects,        /* num_vects for info - adjusted for complex */
                dimension,      /* dimension returned from dloadvector */
                data_type,      /* 1=FLOAT, 2=COMPLEX used as a scale factor */
                max_peaks,      /* max possible peaks = (dimension+1)/2  */
                *ss,            /* slope sign array */
                w1,             /* window to right of "peak" */
                w2,             /* window to left of "peak" */
                w_end,
                flat_start,     /* starting point of a flat region */
                flat_count,     /* count zero slopes found after flat_start */
                max_count,      /* maximum number of peaks found */
                min_count,      /* minimum number of peaks found */
                count,
                i, j, k, l, n, ind,
                val_compare(),
                pos_compare();
  float         **data,         /* vector data returned from dpeakpick */
                sum,
                sumsq,
                temp,
                *sv;           /* slope value array */
  double        w1_avg = 0.0,   /* average slope over window w1 ... 0.0 by jw*/
                w2_avg;         /* average slope over window w2 */


  if(w_size < 1) {
    (void)fprintf(stderr,"ldpeakpick: window (%d) must be greater than 0\n",
                  w_size);
    return(0);
  }
  
  if(hwr < 0) {
    (void)fprintf(stderr,"ldpeakpick: height to width ratio (%f) ",hwr);
    (void)fprintf(stderr,"must be greater than or equal to zero.\n");
    return(0);
  }

 /*
  * load vector array.  
  */
  if(( data = (float **) dload_vector(image, &num_vects, &dimension,
              process_dir)) == NULL) {
    (void)fprintf(stderr,"ldpeakpick: dload_vector failed\n");
    return(0);
  }

 /*
  * If data is complex, dimension returned from dloadvector is 
  * the true dimension.  Set up local variables that double num_vects and
  * half dimension for complex data.
  */
  if(image->data_storage_type == VFF_TYP_COMPLEX) {
    l_vects = num_vects * 2;
    data_type = 2;
  }
  else {
    l_vects = num_vects;
    data_type = 1;
  }

 /*
  * If length of data sequence is less than 3 points, cannot find
  * a peak
  */
  if(dimension < 3) {
    (void)fprintf(stderr,"ldpeakpick: length of sequence must be at ");
    (void)fprintf(stderr,"least 3 points long\n");
    (void)fprintf(stderr,"Sequence length = %d\n", dimension);
    return(0);
  }

  if(!all_flg && (num_peaks < 0) ){
    (void)fprintf(stderr,"ldpeakpick ERROR: number of peaks specified must");
    (void)fprintf(stderr,"be greater zero or greater.\n");
  }

 /*
  * Calculate maximum number of possible peaks.
  */
  max_peaks = (dimension + 1)/2;
  if(!all_flg && (num_peaks > max_peaks) ) {
    (void)fprintf(stderr,"ldpeakpick: specified number of peaks too large\n");
    (void)fprintf(stderr,"     specified:         %d\n",num_peaks);
    (void)fprintf(stderr,"     maximum possible:  %d\n",max_peaks);
    return(0);
  }

 /*
  * allocate memory for peak information array.  
  */

  *peak_info = (struct p_info **) malloc((unsigned int)sizeof(struct p_info *) 
              * l_vects);
  if(*peak_info == NULL)
    (void)fprintf(stderr,"%s: cannot allocate for peak_info array\n",program);

  for(i=0; i<l_vects; i++) {
    (*peak_info)[i] = (struct p_info *) malloc((unsigned int)sizeof 
                      (struct p_info) * max_peaks);
    if((*peak_info)[i] == NULL) 
      (void)fprintf(stderr,"%s: cannot allocate for peak_info array\n",program);
  }

 /*
  * allocate memory for vector information array
  */

  *vect_info = (struct v_info *) malloc((unsigned int)sizeof(struct v_info) 
                * l_vects); 
  if(*vect_info == NULL)
    (void)fprintf(stderr,"%s: cannot allocate for vect_info array\n",program);

 /* 
  * Allocate for temp slope sign and slope value arrays
  */
  if((ss = (int *) malloc((unsigned int)sizeof(int) * (dimension))) == NULL)
    (void)fprintf(stderr,"%s: cannot allocate for slope sign array\n",program);
  if((sv = (float *)malloc((unsigned int)sizeof(float)*(dimension))) == NULL)
    (void)fprintf(stderr,"%s: cannot allocate for slope value array\n",program);

 /*
  * If find troughs is selected, invert data by multiplying by -1
  */
  if(!peak_flg) {
    for(i=0; i<num_vects; i++) {
      for(j=0; j<dimension*data_type; j++) {
        data[i][j] *= -1;
      }
    }
  }


 /*
  * First find all peaks in each vector (handle real and imaginary
  * components of complex data separately).
  */

  switch(image->data_storage_type) {
    case VFF_TYP_FLOAT:
      max_count = 0;
      min_count = max_peaks;
      for(i=0; i<l_vects; i++) {  /* loop through all vectors */
        /* set up slope value and sign arrays for vector i */
        for(j=1; j<dimension; j++) {
          sv[j] = data[i][j] - data[i][j-1];
          if(sv[j] > 0)
            ss[j] = 1;
          else if (sv[j] < 0)
            ss[j] = -1;
          else
            ss[j] = 0;
        }

        /* loop initialization */
        w1=0;
        w2=0;
        count=0;
        n=0;
        flat_start = 0;
        flat_count = 0;

        /* initialize w1 */
        for(l=1; l<=w_size; l++) {
          w1 += ss[l];
          w1_avg += sv[l];
          if(ss[l]==0) flat_count ++;
        }
        w1_avg /= w_size;

        /* is first point in sequence a peak? */
        if(end_flg) {
          if( w1<=(-w_size+flat_count) && fabs(w1_avg)>=hwr ) {
            (*peak_info)[i][count].p_value = data[i][n];
            (*peak_info)[i][count].p_position = n;
            count++;
          }
          else data[i][n] = 0.0;
          flat_count = 0;
        }
        else data[i][n] = 0.0;

        /* detect peaks occuring before enough points to complete w2 */
        for(n=1; n<=w_size; n++) {
          w2 += ss[n];
          w2_avg = 0;
          for(l=1; l<=n; l++)
            w2_avg += sv[l];
          w2_avg /= n;
          w1 += ss[n+w_size] - ss[n];
          w1_avg += (sv[n+w_size] - sv[n])/w_size;
          if(ss[n]==0) flat_count ++;
          if(w2>=n && w1<=-w_size && fabs(w1_avg)>=hwr && fabs(w2_avg)>=hwr) {
            (*peak_info)[i][count].p_value = data[i][n];
            (*peak_info)[i][count].p_position = n;
            count++;
            flat_start = n;
            flat_count = 0;
          }
          else if(w2>=n && fabs(w2_avg)>=hwr) {
            flat_start = n;
            flat_count = 0;
            data[i][n] = 0.0;
          }
          else if(w1<=-w_size && fabs(w1_avg)>=hwr) {
            if( (n-flat_start) == flat_count ) {
              (*peak_info)[i][count].p_value = data[i][n];
              (*peak_info)[i][count].p_position = n;
              count++;
            }
            else data[i][n] = 0.0;
          }
          else if(ss[n]==0 && n!=1) {
            if (ss[n-1] != 0) {
              flat_start = n;
              flat_count = 0;
              data[i][n] = 0.0;
            }
          }
          else data[i][n] = 0.0;
        }

        /* detect peaks where both windows are complete */
        w_end = dimension - w_size;
        for(n=(w_size+1); n<w_end; n++) {
          w1 += ss[n+w_size] - ss[n];
          w1_avg += (sv[n+w_size] - sv[n])/w_size;
          w2 += ss[n] - ss[n-w_size];
          w2_avg += (sv[n] - sv[n-w_size])/w_size;
          if(ss[n]==0) flat_count ++;
          if(w1<=-w_size && w2>=w_size && fabs(w1_avg)>=hwr 
                         && fabs(w2_avg)>=hwr) {
            (*peak_info)[i][count].p_value = data[i][n];
            (*peak_info)[i][count].p_position = n;
            flat_start = n;
            flat_count = 0;
            count++;
          }
          else if(w2>=w_size && fabs(w2_avg)>=hwr) {
            flat_start = n;
            flat_count = 0;
            data[i][n] = 0.0;
          }
          else if(w1<=-w_size && fabs(w1_avg)>=hwr) {
            if( (n-flat_start) == flat_count ) {
              (*peak_info)[i][count].p_value = data[i][n];
              (*peak_info)[i][count].p_position = n;
              count++;
            }
            else data[i][n] = 0.0;
          }
          else if(ss[n]==0 && ss[n-1] != 0) {
            flat_start = n;
            flat_count = 0;
            data[i][n] = 0.0;
          }
          else data[i][n] = 0.0;
        }

        /* detect peaks occuring when not enough points to complete w1 */
        for(n=w_end; n<(dimension-1); n++) {
          w1 -= ss[n];
          w1_avg = 0;
          for (l=n; l<dimension; l++)
            w1_avg += sv[l];
          w1_avg /= dimension - n;
          w2 += ss[n] - ss[n-w_size];
          w2_avg += (sv[n] - sv[n-w_size])/w_size;
          if(ss[n]==0) flat_count ++;
          if(w1<=(n-dimension+1) && w2>=w_size && fabs(w1_avg)>=hwr 
                             && fabs(w2_avg)>=hwr) {
            (*peak_info)[i][count].p_value = data[i][n];
            (*peak_info)[i][count].p_position = n;
            flat_start = n;
            flat_count = 0;
            count++;
          }
          else if(w2>=w_size && fabs(w2_avg)>=hwr) {
            flat_start = n;
            flat_count = 0;
            data[i][n] = 0.0;
          }
          else if(w1<=(n-dimension+1) && fabs(w1_avg)>=hwr) {
            if( (n-flat_start) == flat_count ) {
              (*peak_info)[i][count].p_value = data[i][n];
              (*peak_info)[i][count].p_position = n;
              count++;
            }
            else data[i][n] = 0.0;
          }
          else if(ss[n]==0 && ss[n-1] != 0) {
            flat_start = n;
            flat_count = 0;
            data[i][n] = 0.0;
          }
          else data[i][n] = 0.0;
        }

        /* is last point in sequence a peak? */
        n=dimension - 1;
        if(end_flg) {
          flat_count = 0;
          for(l=(n-w_size); l<n; l++) {
            if(ss[l] == 0) flat_count++;
          }
          w2 += ss[n] - ss[n-w_size];
          w2_avg += (sv[n] - sv[n-w_size])/w_size;
          if( w2>=(w_size-flat_count) && fabs(w2_avg)>=hwr ) {
            (*peak_info)[i][count].p_value = data[i][n];
            (*peak_info)[i][count].p_position = n;
            count++;
          }
          else data[i][n] = 0.0;
        }
        else data[i][n] = 0.0;

        if(count > max_count)
          max_count = count;
        if(count < min_count)
          min_count = count;
        (*vect_info)[i].p_count = count;

      } /* end of loop through vectors */

      /* If find maximum amplitude peaks was specified, have to sort the 
         array by maximum value to select the "num_peaks" maximum peaks. */
      if(max_flg && !all_flg)
      {
        for(i=0; i<l_vects; i++) {
          (void)qsort((char *)(char *)(*peak_info)[i], (*vect_info)[i].p_count, 
                        sizeof(struct p_info), val_compare);
        }
      }
  
      /* zero out peak values past the first "num_peaks" in list" */
      if(!all_flg) {
        for(i=0; i<l_vects; i++) {
          for(j=num_peaks; j<(*vect_info)[i].p_count; j++) {
            ind = (*peak_info)[i][j].p_position;
            data[i][ind] = 0.0;
          }
          if((*vect_info)[i].p_count > num_peaks)
            (*vect_info)[i].p_count = num_peaks;
        }
      }

      /* If values were sorted above, need to resort back to positional order */
      if(max_flg && !all_flg) {
        for(i=0; i<l_vects; i++) {
          (void)qsort((char *)(char *)(*peak_info)[i], (*vect_info)[i].p_count, 
                        sizeof(struct p_info), pos_compare);
        }
      }

      /* calculate mean and variance of peak to peak distances */ 
      for(i=0; i<l_vects; i++) {
        if((*vect_info)[i].p_count > 1) {
          sum = 0;
          sumsq=0;
          for(j=0; j<((*vect_info)[i].p_count - 1); j++) {
            temp = (*peak_info)[i][j+1].p_position - 
                          (*peak_info)[i][j].p_position;
            sum += temp;
            temp *= temp;
            sumsq += temp;
          }
          temp = sum/((*vect_info)[i].p_count - 1);
          (*vect_info)[i].p2p_dist_mean = temp;
          temp *= temp;
          (*vect_info)[i].p2p_dist_var = (sumsq/((*vect_info)[i].p_count - 1)) 
                                                - temp;
        }
        else {
          (*vect_info)[i].p2p_dist_mean = 0.0;
          (*vect_info)[i].p2p_dist_var  = 0.0;
        }
      }
          
    break;

    case VFF_TYP_COMPLEX:
      max_count = 0;
      min_count = max_peaks;

      for(i=0; i<num_vects; i++) {  /* loop through all vectors */
        /* set up slope value and sign arrays for vector i */
        for(k=0; k<2; k++) { /* first do all reals then all imaginaries */
          ind = 1;
          for(j=(k+2); j<dimension*data_type; j+=2) {
            sv[ind] = data[i][j] - data[i][j-2];
            if(sv[ind] > 0)
              ss[ind] = 1;
            else if (sv[ind] < 0)
              ss[ind] = -1;
            else
              ss[ind] = 0;
            ind++;
          }

          /* loop initialization */
          w1=0;
          w2=0;
          n=0;
          count=0;
          flat_start = 0;
          flat_count = 0;
          ind = 2*i + k;

          /* initialize w1 */
          for(l=1; l<=w_size; l++) {
            w1 += ss[l];
            w1_avg += sv[l];
            if(ss[l]==0) flat_count ++;
          }
          w1_avg /= w_size;

          /* is first point in sequence a peak? */
          if(end_flg) {
            if( w1<=(-w_size+flat_count) && fabs(w1_avg)>=hwr ) {
              (*peak_info)[ind][count].p_value = data[i][(2*n)+k];
              (*peak_info)[ind][count].p_position = n;
              count++;
            }
            else data[i][(2*n)+k] = 0.0;
            flat_count = 0;
          }
          else data[i][(2*n)+k] = 0.0;


          /* detect peaks occuring before enough points to complete w2 */
          for(n=1; n<=w_size; n++) {
            w2 += ss[n];
            w2_avg = 0;
            for(l=1; l<=n; l++)
              w2_avg += sv[l];
            w2_avg /= n;
            w1 += ss[n+w_size] - ss[n];
            w1_avg += (sv[n+w_size] - sv[n])/w_size;
            if(ss[n]==0) flat_count ++;
            if(w2>=n && w1<=-w_size && fabs(w1_avg)>=hwr && fabs(w2_avg)>=hwr) {
              (*peak_info)[ind][count].p_value = data[i][(2*n)+k];
              (*peak_info)[ind][count].p_position = n;
              count++;
              flat_start = n;
              flat_count = 0;
            }
            else if(w2>=n && fabs(w2_avg)>=hwr) {
              flat_start = n;
              flat_count = 0;
              data[i][(2*n)+k] = 0.0;
            }
            else if(w1<=-w_size && fabs(w1_avg)>=hwr) {
              if( (n-flat_start) == flat_count ) {
                (*peak_info)[ind][count].p_value = data[i][(2*n)+k];
                (*peak_info)[ind][count].p_position = n;
                count++;
              }
              else data[i][(2*n)+k] = 0.0;
            }
            else if(ss[n]==0 && n!=1) {
              if (ss[n-1] != 0) {
                flat_start = n;
                flat_count = 0;
                data[i][(2*n)+k] = 0.0;
              }
            }
            else data[i][(2*n)+k] = 0.0;
          }

          /* detect peaks where both windows are complete */
          w_end = dimension - w_size;
          for(n=(w_size+1); n<w_end; n++) {
            w1 += ss[n+w_size] - ss[n];
            w1_avg += (sv[n+w_size] - sv[n])/w_size;
            w2 += ss[n] - ss[n-w_size];
            w2_avg += (sv[n] - sv[n-w_size])/w_size;
            if(ss[n]==0) flat_count ++;
            if(w1<=-w_size && w2>=w_size && fabs(w1_avg)>=hwr 
                           && fabs(w2_avg)>=hwr) {
              (*peak_info)[ind][count].p_value = data[i][(n*2)+k];
              (*peak_info)[ind][count].p_position = n;
              flat_start = n;
              flat_count = 0;
              count++;
            }
            else if(w2>=w_size && fabs(w2_avg)>=hwr) {
              flat_start = n;
              flat_count = 0;
              data[i][(n*2)+k] = 0.0;
            }
            else if(w1<=-w_size && fabs(w1_avg)>=hwr) {
              if( (n-flat_start) == flat_count ) {
                (*peak_info)[ind][count].p_value = data[i][(n*2)+k];
                (*peak_info)[ind][count].p_position = n;
                count++;
              }
              else data[i][(n*2)+k] = 0.0;
            }
            else if(ss[n]==0 && ss[n-1] != 0) {
              flat_start = n;
              flat_count = 0;
              data[i][(n*2)+k] = 0.0;
            }
            else data[i][(n*2)+k] = 0.0;
          }

          /* detect peaks occuring when not enough points to complete w1 */
          for(n=w_end; n<(dimension-1); n++) {
            w1 -= ss[n];
            w1_avg = 0;
            for (l=n; l<dimension; l++)
              w1_avg += sv[l];
            w1_avg /= dimension - n;
            w2 += ss[n] - ss[n-w_size];
            w2_avg += (sv[n] - sv[n-w_size])/w_size;
            if(ss[n]==0) flat_count ++;
            if(w1<=(n-dimension+1) && w2>=w_size && fabs(w1_avg)>=hwr 
                               && fabs(w2_avg)>=hwr) {
              (*peak_info)[ind][count].p_value = data[i][(2*n)+k];
              (*peak_info)[ind][count].p_position = n;
              flat_start = n;
              flat_count = 0;
              count++;
            }
            else if(w2>=w_size && fabs(w2_avg)>=hwr) {
              flat_start = n;
              flat_count = 0;
              data[i][(2*n)+k] = 0.0;
            }
            else if(w1<=(n-dimension+1) && fabs(w1_avg)>=hwr) {
              if( (n-flat_start) == flat_count ) {
                (*peak_info)[ind][count].p_value = data[i][(2*n)+k];
                (*peak_info)[ind][count].p_position = n;
                count++;
              }
              else data[i][(2*n)+k] = 0.0;
            }
            else if(ss[n]==0 && ss[n-1] != 0) {
              flat_start = n;
              flat_count = 0;
              data[i][(2*n)+k] = 0.0;
            }
            else data[i][(2*n)+k] = 0.0;
          }

          /* is last point in sequence a peak? */
          n=dimension - 1;
          if(end_flg) {
            flat_count = 0;
            for(l=(n-w_size); l<n; l++) {
              if(ss[l] == 0) flat_count++;
            }
            w2 += ss[n] - ss[n-w_size];
            w2_avg += (sv[n] - sv[n-w_size])/w_size;
            if( w2>=(w_size-flat_count) && fabs(w2_avg)>=hwr ) {
              (*peak_info)[ind][count].p_value = data[i][(2*n)+k];
              (*peak_info)[ind][count].p_position = n;
              count++;
            }
            else data[i][(2*n)+k] = 0.0;
          }
          else data[i][(2*n)+k] = 0.0;

          if(count > max_count)
            max_count = count;
          if(count < min_count)
            min_count = count;
          (*vect_info)[ind].p_count = count;

        } /* end of k loop */

      } /* end of loop through vectors */

      /* If find maximum amplitude peaks was specified, have to sort the 
         array by maximum value to select the "num_peaks" maximum peaks. */
      if(max_flg && !all_flg)
      {
        for(i=0; i<l_vects; i++) {
          (void)qsort((char *)(*peak_info)[i], (*vect_info)[i].p_count, 
                        sizeof(struct p_info), val_compare);
        }
      }
  
      /* zero out peak values past the first "num_peaks" in list" */
      if(!all_flg) {
        k=1;
        for(i=0; i<l_vects; i++) {
          for(j=num_peaks; j<(*vect_info)[i].p_count; j++) { 
            if(k>0)  /* real component */
              ind = 2*(*peak_info)[i][j].p_position;
            else   /* imaginary component */
              ind = 2*(*peak_info)[i][j].p_position + 1;
            data[(int)(i/2)][ind] = 0.0;
          }
          if((*vect_info)[i].p_count > num_peaks)
            (*vect_info)[i].p_count = num_peaks;
          k *= -1;
        }
      }

      /* If values were sorted above, need to resort back to positional order */
      if(max_flg && !all_flg) {
        for(i=0; i<l_vects; i++) {
          (void)qsort((char *)(*peak_info)[i], (*vect_info)[i].p_count, 
                        sizeof(struct p_info), pos_compare);
        }
      }

      /* calculate mean and variance of peak to peak distances */ 
      for(i=0; i<l_vects; i++) {
        if((*vect_info)[i].p_count > 1) {
          sum = 0;
          sumsq=0;
          for(j=0; j<((*vect_info)[i].p_count - 1); j++) {
            temp = (*peak_info)[i][j+1].p_position - 
                          (*peak_info)[i][j].p_position;
            sum += temp;
            temp *= temp;
            sumsq += temp;
          }
          temp = sum/((*vect_info)[i].p_count - 1);
          (*vect_info)[i].p2p_dist_mean = temp;
          temp *= temp;
          (*vect_info)[i].p2p_dist_var = (sumsq/((*vect_info)[i].p_count - 1)) 
                                                - temp;
        }
        else {
          (*vect_info)[i].p2p_dist_mean = 0.0;
          (*vect_info)[i].p2p_dist_var  = 0.0;
        }
      }
          
    break;

    default:
      (void)fprintf(stderr,"ldpeakpick: unknown data type\n");
      return(0);
    
  } /* end of data_storage_type case statement for finding peaks */

    
  if(!all_flg && (min_count < num_peaks)) {
    (void)fprintf(stderr,"ldpeakpick:\n");
    (void)fprintf(stderr,"WARNING: one or more vectors contains less");
    (void)fprintf(stderr,"than specified number of peaks\n");
    (void)fprintf(stderr,"     specified number of peaks = %d\n",num_peaks);
    (void)fprintf(stderr,"     minimum number found      = %d\n",min_count);
  }

  /* reallocate peak_info array */
  if(all_flg) max_peaks = max_count;
  else max_peaks = num_peaks;
  for(i=0; i<l_vects; i++) {
    (*peak_info)[i] = (struct p_info *) realloc((*peak_info)[i], 
                 (unsigned int)sizeof(struct p_info) * max_peaks);
    if((*peak_info)[i] == NULL) 
      (void)fprintf(stderr,"%s: cannot reallocate peak_info array\n",program);
  }
  

 /*
  * If find troughs is selected, invert peak values by multiplying by -1
  */
  if (!peak_flg) {
    for (i=0; i<num_vects; i++) {
      for (j=0; j<dimension*data_type; j++)
        data[i][j] *= -1;
    }
    for (i=0; i<l_vects; i++) {
      for (j=0; j<(*vect_info)[i].p_count; j++)
        (*peak_info)[i][j].p_value *= -1;
    }
  }

 /*
  * Unload data back into image format (dunload_vector)
  */

  if (!dunload_vector((char **)data, image, image->data_storage_type,
                      num_vects, dimension, process_dir) ) {
    (void)fprintf(stderr,"ldpeakpick: dunload_vector failed\n");
    return(0);
  }

  for (i=0; i<num_vects; i++)
    free((char*)data[i]);
  free((char**)data);

  return(1);

} /* end of ldpeakpick */


/*
 * val_compare sorts the structure by value.  First element in resulting 
 * array will be largest value, and last element will be smallest value.
 */
int val_compare(arg1, arg2)
struct p_info *arg1, *arg2;
{
  if(arg1->p_value < arg2->p_value)
    return(1);
  if(arg1->p_value > arg2->p_value)
    return(-1);
    return(0);
}

/*
 * pos_compare sorts the structure by position.  First element in resulting 
 * array will correspond to smallest position value, and last element will 
 * correspond to largest position value.
 */
int pos_compare(arg1, arg2)
struct p_info *arg1, *arg2;
{
  if(arg1->p_position < arg2->p_position)
    return(-1);
  if(arg1->p_position > arg2->p_position)
    return(1);
    return(0);
}


/* -library_code_end */
