package yi.math;

public class Convolution {
  public final static int SUM_TO_ZERO = 0;
  public final static int BALANCE_WEIGHT = 1;
  public final static int REFLECTION = 2;
  public int BOUNDARY_PROCESS_METHOD = SUM_TO_ZERO;

  public static double[][] convolution(
                             double[][] in, int in_w, int in_h,
                             double[][] filter, int filter_w, int filter_h,
                             int boundary_process_method) {
    double[][] out = new double[in_h][in_w];
    if (boundary_process_method == SUM_TO_ZERO) {
      for (int y=0; y<in_h; y++) {
        for (int x=0; x<in_w; x++) {
          // scan
          int y_in = y - filter_h / 2;
          double sum = 0, filtered_sum = 0, filter_sum = 0, filter_num = 0;
          for (int sy=0; sy<filter_h; sy++, y_in++) {
            if (y_in<0||y_in>=in_h) continue;
            int x_in = x - filter_w/2;
            for (int sx=0; sx<filter_w; sx++,x_in++) {
              if (x_in<0||x_in>=in_w) continue;
              sum += in[y_in][x_in];
              filter_sum += filter[sy][sx];
              filter_num ++;
              filtered_sum += in[y_in][x_in] * filter[sy][sx];
            }
          }
          out[y][x] = filtered_sum - sum * filter_sum / filter_num;
        } // x
      } // y
      return out;
    }
    else if (boundary_process_method == BALANCE_WEIGHT) {
      for (int y=0; y<in_h; y++) {
        for (int x=0; x<in_w; x++) {
          // scan
          int y_in = y - filter_h / 2;
          double sum = 0, filtered_sum = 0, filter_sum = 0, filter_num = 0;
          for (int sy=0; sy<filter_h; sy++, y_in++) {
            if (y_in<0||y_in>=in_h) continue;
            int x_in = x - filter_w/2;
            for (int sx=0; sx<filter_w; sx++,x_in++) {
              if (x_in<0||x_in>=in_w) continue;
              sum += in[y_in][x_in];
              filter_sum += filter[sy][sx];
              filter_num ++;
              filtered_sum += in[y_in][x_in] * filter[sy][sx];
            }
          }
          out[y][x] = filtered_sum - sum * filter_sum / filter_num;
          out[y][x] *= filter_w * filter_h / (double) filter_num;
        } // x
      } // y
      return out;
    }
    else if (boundary_process_method == REFLECTION) {
      for (int y=0; y<in_h; y++) {
        for (int x=0; x<in_w; x++) {
          // scan
          int y_in = y - filter_h / 2;
          double sum = 0, filtered_sum = 0, filter_sum = 0, filter_num = 0;
          for (int sy=0; sy<filter_h; sy++, y_in++) {
            if (y_in<0) y_in = -y_in;
            if (y_in>=in_h) y_in = in_h - 1 - (y_in - in_h);
            int x_in = x - filter_w/2;
            for (int sx=0; sx<filter_w; sx++,x_in++) {
              if (x_in<0) x_in = -x_in;
              if (x_in>=in_w) x_in = in_w - 1 - (x_in-in_w);
              sum += in[y_in][x_in];
            }
          }
        } // x
      } // y
      return out;
    }
    else return null;
  }

}