package nxm.dsp.lib;

/**
  Implements morphological filter elements.

  @author Jeff Schoen
*/
public class Morph {

  /** This constructor was automatically provided by Java,
   *  but it is unneeded as all the methods are static methods.
   *  @deprecated now that Javadocs are required, it is temporarily
   *              being explicitely provided but deprecated
   */
  @Deprecated
  public Morph () {}

  /** Minimum
   *  @param vi   input data buffer
   *  @param vo   output data buffer
   *  @param ntot number of elements in input data buffer
   *  @param nw   window length
   */
  public static void min (double[] vi, double[] vo, int ntot, int nw) {
    int i, il, ir, nl=nw/2, nr=nw-nl;
    // section A
    for (i=il=0,ir=nr; i<nl; i++,ir++) vo[i] = min(vi,0,ir);
    // section C
    for (il--,ir--; i<ntot-nr; i++,il++,ir++) {
      if (vi[ir] < vo[i-1]) vo[i] = vi[ir];                 // next < prev MIN
      else if (il>=0 && vi[il] > vo[i-1]) vo[i] = vo[i-1];  // last > prev MIN
      else vo[i] = min(vi,il+1,nw);                         // recompute
    }
    // section E
    for (il++,ir++; i<ntot; i++,il++) vo[i] = min(vi,il,ir-il);
  }

  private static double min (double[] v, int i, int n) {
    double w = v[i];
    n+=i; for (i++; i<n; i++) { if (v[i]<w) w=v[i]; }
    return w;
  }

  /** Maximum Morphological operator
   *  @param vi   input data buffer
   *  @param vo   output data buffer
   *  @param ntot number of elements in input data buffer
   *  @param nw   window length
   */
  public static void max (double[] vi, double[] vo, int ntot, int nw) {
    int i, il, ir, nl=nw/2, nr=nw-nl;
    // section A
    for (i=il=0,ir=nr; i<nl; i++,ir++) vo[i] = max(vi,0,ir);
    // section C
    for (il--,ir--; i<ntot-nr; i++,il++,ir++) {
      if (vi[ir] > vo[i-1]) vo[i] = vi[ir];                 // next > prev MAX
      else if (il>=0 && vi[il] < vo[i-1]) vo[i] = vo[i-1];  // last < prev MAX
      else vo[i] = max(vi,il+1,nw);                         // recompute
    }
    // section E
    for (il++,ir++; i<ntot; i++,il++) vo[i] = max(vi,il,ir-il);
  }

  private static double max (double[] v, int i, int n) {
    double w = v[i];
    n+=i; for (i++; i<n; i++) { if (v[i]>w) w=v[i]; }
    return w;
  }

  /** Mean Morphological operator
   *  @param vi   input data buffer
   *  @param vo   output data buffer
   *  @param ntot number of elements in input data buffer
   *  @param nw   window length
   */
  public static void mean (double[] vi, double[] vo, int ntot, int nw) {
    int i, il, ir, iw, nl=nw/2, nr=nw-nl;
    double sum=0, scale=1.0/nw;
    // section A
    for (i=il=0,ir=nr; i<nl; i++,ir++) vo[i] = sum(vi,0,ir)/ir;
    // section C
    for (il--,ir--,iw=nw; i<ntot-nr; i++,il++,ir++,iw++) {
      if (iw>=nw) { sum = sum(vi,i-nl,nw); iw=0; }
      else sum = sum - vi[il] + vi[ir];
      vo[i] = sum*scale;
    }
    // section E
    for (il++,ir++; i<ntot; i++,il++) vo[i] = sum(vi,il,ir-il)/(ir-il);
  }

  private static double sum (double[] v, int i, int n) {
    double w = v[i];
    n+=i; for (i++; i<n; i++) w+=v[i];
    return w;
  }

  /** Median Morphological operator
   *  @param vi   input data buffer
   *  @param vo   output data buffer
   *  @param ntot number of elements in input data buffer
   *  @param nw   window length
   */
  public static void median (double[] vi, double[] vo, int ntot, int nw) {
    int i, il, ir, iw, nl=nw/2, nr=nw-nl;
    double median=0; double[] sort = new double[nw];
    // section A
    for (i=il=0,ir=nr; i<nl; i++,ir++) vo[i] = median(vi,0,ir,sort);
    // section C
    for (il--,ir--; i<ntot-nr; i++,il++,ir++) {
//      vo[i] = median(vo[i-1],vi[il],vi[ir],nw,sort,nl);
      vo[i] = median(vi,i-nl,nw,sort);
    }
    // section E
    for (il++,ir++; i<ntot; i++,il++) vo[i] = median(vi,il,ir-il,sort);
  }

  private static double median (double[] v, int i1, int n, double[] sort) {
    int i,j,k;
    for (i=0; i<n; i++) {
      double s = v[i1+i];
      for (j=0; j<i && s>sort[j]; j++);
      for (k=i; k>j; k--) sort[k]=sort[k-1];
      sort[j]=s;
    }
    return sort[n/2];
  }

}
