package yi.iu.processor.texture.gabor;

import yi.util.io.*;
import yi.util.io.input.*;
import yi.util.io.output.*;
import yi.util.config.*;
import yi.util.log.*;


import yi.math.DoubleMatrix;
import java.awt.*;
import yi.iu.processor.color.quan.*;
import yi.iu.processor.channel.*;
import yi.iu.processor.resample.*;
import yi.util.data.*;
import yi.iu.processor.filter.*;
import java.util.*;

/**
 * Class Gabor is for Gabor filtering. Details are given in the paper
 * by
 *
 * Manjunath and Ma, "Texture features for browsing and retrieval of image
 * data", IEEE PAMI, vol 18, no. 8, pp. 837-842, Aug. 1996.
 *
 * Original code was implemented by Ma in C.
 *
 * @author Selim Aksoy
 * @version 1.0 (3/26/2001)
 */
public class Gabor implements Configurable, RstReusable{

  public int    REUSE  = DIU.INPUT_NONE;
  public int    OUTPUT = DOU.OUTPUT_NONE;

  public int    LOG_LEVEL = 3;

  public String[] getParameterNames() {
    return new String[] {
      "LOG_LEVEL",
      "REUSE",
      "OUTPUT"
    };
  }

  static double freqLower = 0.15;	//lower frequency
  static double freqUpper = 0.4;  //upper frequency
  public static int numScales = 3;       //number of scales in generating filters
  public static int numOrients = 4;      //number of orientations in generating filters
  static int removeDC = 1;        //flag for removing the DC component
  public static int maskSize = 10;       //filter mask = 2*maskSize+1 x 2*maskSize+1

  //constructor
  public Gabor() {
    Configuration.configParameters(this);
  }

  public static DoubleMatrix getFiltered(DoubleMatrix image, int numScales, int numOrients ) {

          int maskSizeTotal = 2 * maskSize + 1;

    //allocate space for filtering results
    DoubleMatrix gaborMaskReal = new DoubleMatrix( maskSizeTotal, maskSizeTotal );
    DoubleMatrix gaborMaskImag = new DoubleMatrix( maskSizeTotal, maskSizeTotal );

    DoubleMatrix features = new DoubleMatrix( image.height * image.width,
            numScales * numOrients );

//          System.out.println( "numScales: " + numScales + ", numOrients: " + numOrients );

    //for all scales and orientations
    int filter_index = 0;
    for(int s = 0; s < numScales; s++ ) {
      for(int n = 0; n < numOrients; n++, filter_index++ ) {

        System.out.println( "s: " + s + ", n: " + n );

        //form the Gabor filter
        gaborFilter( gaborMaskReal, gaborMaskImag, numScales, numOrients, s + 1, n + 1 );

        double real_weight = yi.stat.StatUtil.absoluteSum(gaborMaskReal.data);
        double imag_weight = yi.stat.StatUtil.absoluteSum(gaborMaskImag.data);

        double[][] real = yi.math.Convolution.convolution(
                            image.data, image.width, image.height,
                            gaborMaskReal.data, maskSizeTotal, maskSizeTotal,
                            yi.math.Convolution.SUM_TO_ZERO);

        double[][] imag = yi.math.Convolution.convolution(
                            image.data, image.width, image.height,
                            gaborMaskImag.data, maskSizeTotal, maskSizeTotal,
                            yi.math.Convolution.SUM_TO_ZERO);

        //compute the magnitude of the complex
        //filtering results
        int data_index = 0;
        for(int r = 0; r < image.height; r++ ) {
                for(int c = 0; c < image.width; c++, data_index++ ) {
                        features.data[ data_index ][ filter_index ] = Math.sqrt(
                                Math.pow( real[ r ][ c ]/real_weight, 2.0 ) +
                                Math.pow( imag[ r ][ c ]/imag_weight, 2.0 ) );
                }
        }

      }
    }
    return features;

  }

  //main method to filter the input image using multiple Gabor filters
  public static DoubleMatrix getProcessed( DoubleMatrix image, int numScales, int numOrients ) {

          int r, c, s, n, i, j;	//array indices
          int height, width;		//original image size
          int rowSize, colSize;	//filtered image size
          int maskSizeTotal = 2 * maskSize + 1;

          //find subsampled image size
          rowSize = (int) Math.pow( 2.0, Math.ceil(
                  Math.log( (double) image.height ) / Math.log( 2.0 ) ) );
          colSize = (int) Math.pow( 2.0, Math.ceil(
                  Math.log( (double) image.width ) / Math.log( 2.0 ) ) );
          //System.out.println( "rowSize: " + rowSize + " colSize: " + colSize );

          //allocate space for FFT results
          DoubleMatrix imageReal = new DoubleMatrix( rowSize, colSize );
          DoubleMatrix imageImag = new DoubleMatrix( rowSize, colSize );
          DoubleMatrix fftReal = new DoubleMatrix( rowSize, colSize );
          DoubleMatrix fftImag = new DoubleMatrix( rowSize, colSize );

          //compute FFT of the input image
          for( r = 0; r < image.height; r++ )
                  for( c = 0; c < image.width; c++ )
                          imageReal.data[ r ][ c ] = image.data[ r ][ c ];
          DoubleMatrix.fft2( fftReal, fftImag, imageReal, imageImag );

          //allocate space for filtering results
          DoubleMatrix gaborMaskReal = new DoubleMatrix( maskSizeTotal, maskSizeTotal );
          DoubleMatrix gaborMaskImag = new DoubleMatrix( maskSizeTotal, maskSizeTotal );
          DoubleMatrix gaborReal = new DoubleMatrix( rowSize, colSize );
          DoubleMatrix gaborImag = new DoubleMatrix( rowSize, colSize );
          DoubleMatrix features = new DoubleMatrix( image.height * image.width,
                  numScales * numOrients );

//                System.out.println( "numScales: " + numScales + ", numOrients: " + numOrients );

          //for all scales and orientations
          for( s = 0; s < numScales; s++ ) {
                  for( n = 0; n < numOrients; n++ ) {

//                                System.out.println( "s: " + s + ", n: " + n );

                          //form the Gabor filter
                          gaborFilter( gaborMaskReal, gaborMaskImag, numScales, numOrients, s + 1, n + 1 );
//                                System.out.println("max,min,sum of REAL= " +
//                                  yi.stat.Util.maximum(gaborMaskReal.data) + "," +
//                                  yi.stat.Util.minimum(gaborMaskReal.data) + "," +
//                                  yi.stat.Util.sum(gaborMaskReal.data));
//                                System.out.println("max,min,sum of IMAG= " +
//                                  yi.stat.Util.maximum(gaborMaskImag.data) + "," +
//                                  yi.stat.Util.minimum(gaborMaskImag.data) + "," +
//                                  yi.stat.Util.sum(gaborMaskImag.data));
//                                yi.util.io.OutputResult.print("gaborMaskReal" + s + "." + n, gaborMaskReal.data, maskSizeTotal, maskSizeTotal);
//                                yi.util.io.OutputResult.print("gaborMaskImag" + s + "." + n, gaborMaskImag.data, maskSizeTotal, maskSizeTotal);
                          //compute FFT of Gabor filter
        double real_weight = yi.stat.StatUtil.absoluteSum(gaborMaskReal.data);
        double imag_weight = yi.stat.StatUtil.absoluteSum(gaborMaskImag.data);
//        System.out.println("real weight= " + real_weight);
//        System.out.println("imag weight= " + imag_weight);
                          DoubleMatrix.fft2( gaborReal, gaborImag,
                                  DoubleMatrix.copy( gaborMaskReal, rowSize, colSize,
                                          0, 0, 0, 0, 2 * maskSize, 2 * maskSize ),
                                  DoubleMatrix.copy( gaborMaskImag, rowSize, colSize,
                                          0, 0, 0, 0, 2 * maskSize, 2 * maskSize ) );

                          //convolution is done as multiplication
                          //in the frequency domain
                          imageReal = DoubleMatrix.subtract(
                                  DoubleMatrix.multiply( gaborReal, fftReal ),
                                  DoubleMatrix.multiply( gaborImag, fftImag ) );
                          imageImag = DoubleMatrix.add(
                                  DoubleMatrix.multiply( gaborReal, fftImag ),
                                  DoubleMatrix.multiply( gaborImag, fftReal ) );

                          //take inverse FFT and get the
                          //convolution back in time domain
                          DoubleMatrix.ifft2( gaborReal, gaborImag, imageReal, imageImag );
                          imageReal = DoubleMatrix.shift( gaborReal, maskSize );
                          imageImag = DoubleMatrix.shift( gaborImag, maskSize );

                          //compute the magnitude of the complex
                          //filtering results
                          i = 0;
                          j = s * numOrients + n;
                          for( r = 0; r < image.height; r++ ) {
                                  for( c = 0; c < image.width; c++ ) {
                                          features.data[ i ][ j ] = Math.sqrt(
                                                  Math.pow( imageReal.data[ r ][ c ]/real_weight, 2.0 ) +
                                                  Math.pow( imageImag.data[ r ][ c ]/imag_weight, 2.0 ) );
                                          i++;
                                  }
                          }

//                          // smooth first
//                          double[] to_smooth = new double[image.height * image.width];
//                          i=0;
//                          for( r = 0; r < image.height; r++ ) {
//                                  for( c = 0; c < image.width; c++, i++ ) {
//                                          to_smooth[ i ] = Math.sqrt(
//                                                  Math.pow( imageReal.data[ r ][ c ]/real_weight, 2.0 ) +
//                                                  Math.pow( imageImag.data[ r ][ c ]/imag_weight, 2.0 ) );
//                                  }
//                          }
//                          GaussianSmooth smoother = new GaussianSmooth(2);
//                          double[] smoothed = smoother.getProcessed(to_smooth, image.width,image.height);
//
//                          // assign the smoothed to features
//                          j = s * numOrients + n;
//                          for (i=0; i<smoothed.length; i++) features.data[i][j] = smoothed[i];

                  }
          }

          return features;

  }

  //method to generate Gabor filter masks in time domain
  public static void gaborFilter( DoubleMatrix gaborReal, DoubleMatrix gaborImag,  int numScales, int numOrients, int s, int n ) {

          double freqBase, a, u0, z, v, x0, y0, g, t1, t2, m;
          int x, y, side;

          freqBase = freqUpper / freqLower;
          a = Math.pow( freqBase, 1.0 / (double) ( numScales - 1 ) );

          u0 = freqUpper / Math.pow( a, (double) ( numScales - s ) );

          v = Math.pow( 0.6 / freqUpper *
                  Math.pow( a, (double) ( numScales - s ) ), 2.0 );

          t1 = Math.cos( (double) Math.PI / numOrients * ( n - 1.0 ) );
          t2 = Math.sin( (double) Math.PI / numOrients * ( n - 1.0 ) );

          side = (int) ( gaborReal.height - 1 ) / 2;

          for( x = 0; x < ( 2 * side + 1 ); x++ ) {
                  for( y = 0; y < ( 2 * side + 1 ); y++ ) {
                          x0 = (double) ( x - side) * t1 + (double) ( y - side ) * t2;
                          y0 = (double) - ( x - side ) * t2 + (double) ( y - side ) * t1;
                          g = 1.0 / ( 2.0 * Math.PI * v ) *
                                  Math.pow( a, (double) ( numScales - s ) ) *
                                  Math.exp( - 0.5 * ( x0 * x0 + y0 * y0 ) / v );
                          gaborReal.data[x][y] = g * Math.cos( 2.0 * Math.PI * u0 * x0 );
                          gaborImag.data[x][y] = g * Math.sin( 2.0 * Math.PI * u0 * x0 );
                  }
          }

          //if removeDC flag is set, then remove the average value
          //from the real part of Gabor
          if( removeDC == 1 ) {

                  m = 0;
                  for( x = 0; x < ( 2 * side + 1 ); x++ )
                          for( y = 0; y < ( 2 * side + 1 ); y++ )
                                  m += gaborReal.data[x][y];

                  m /= Math.pow( (double) 2.0 * side + 1, 2.0 );

                  for( x = 0; x < ( 2 * side + 1 ); x++ )
                          for( y = 0; y < ( 2 * side + 1 ); y++ )
                                  gaborReal.data[x][y] -= m;

          }

  }


  public Object reuse(int reuse_level) {
    if (reuse_level != DIU.INPUT_NONE) {
      try {
        // only serializer is supported for this class
        if (reuse_level != DIU.INPUT_SER) {
          throw new DataInputException("reuse method " + reuse_level +
            "is not supported by Gabor");
        }
        DoubleMatrix rst = (DoubleMatrix)DIU.in(
            reuse_level,
            DOU.shortName(this.getClass().getName()),
            "GaborRst",
            null);
        return rst;
      }
      catch (Exception e) {
        LogHandler.log(this, LOG_LEVEL, 1, "Can't Reuse Energy Measures:\n" + e);
      }
    }
    return null;
  }

  void output(int output_level, DoubleMatrix ds){
    if (output_level != DOU.OUTPUT_NONE) {
      try {
        // only serializer is supported for this class
        if (output_level != DOU.OUTPUT_SER) {
          throw new DataOutputException("output method " + output_level +
            "is not supported by Gabor");
        }
        DOU.out(
            output_level,
            DOU.shortName(this.getClass().getName()),
            "GaborRst",
            ds);
      } catch (DataOutputException e) {
        LogHandler.log(this, LOG_LEVEL, 1, e.toString());
      }
    }
  }

  public DoubleMatrix getReusableGaborRst(DoubleMatrix image, int numScales, int numOrients ) {
      DoubleMatrix rst = (DoubleMatrix)reuse(REUSE);

      if (rst == null) {
        rst = Gabor.getProcessed(image, numScales, numOrients);
      }

      output(OUTPUT, rst);

      return rst;
  }

  public static class Filters {
    /**
     *
     */
    public static void main(String[] args) {

      int r, c, s, n, i, j;	//array indices
      int height, width;		//original image size
      int rowSize, colSize;	//filtered image size
      int maskSizeTotal = 2 * maskSize + 1;

      //allocate space for filtering results
      DoubleMatrix gaborMaskReal = new DoubleMatrix( maskSizeTotal, maskSizeTotal );
      DoubleMatrix gaborMaskImag = new DoubleMatrix( maskSizeTotal, maskSizeTotal );

      System.out.println( "numScales: " + numScales + ", numOrients: " + numOrients );

      //for all scales and orientations
      for( s = 0; s < numScales; s++ ) {
          for( n = 0; n < numOrients; n++ ) {

              System.out.println( "s: " + s + ", n: " + n );

              //form the Gabor filter
              gaborFilter( gaborMaskReal, gaborMaskImag, numScales, numOrients, s + 1, n + 1 );
//                                System.out.println("max,min,sum of REAL= " +
//                                  yi.stat.Util.maximum(gaborMaskReal.data) + "," +
//                                  yi.stat.Util.minimum(gaborMaskReal.data) + "," +
//                                  yi.stat.Util.sum(gaborMaskReal.data));
//                                System.out.println("max,min,sum of IMAG= " +
//                                  yi.stat.Util.maximum(gaborMaskImag.data) + "," +
//                                  yi.stat.Util.minimum(gaborMaskImag.data) + "," +
//                                  yi.stat.Util.sum(gaborMaskImag.data));
              double real_weight = yi.stat.StatUtil.absoluteSum(gaborMaskReal.data);
              double imag_weight = yi.stat.StatUtil.absoluteSum(gaborMaskImag.data);
              System.out.println("real weight= " + real_weight);
              System.out.println("imag weight= " + imag_weight);

              yi.util.io.OutputResult.print("gaborMaskReal" + s + "." + n, gaborMaskReal.data, maskSizeTotal, maskSizeTotal);
              yi.util.io.OutputResult.print("gaborMaskImag" + s + "." + n, gaborMaskImag.data, maskSizeTotal, maskSizeTotal);
          }
      }
    }
  }

  /**
   * Feed Gabor feature extraction processor with luminance data
   * @see #GABOR_DATA_FEED
   */
  public final static int GABOR_USE_LUMINANCE = 0;
  /**
   * Feed Gabor feature extraction processor with average of H,S,I.  This method
   * is used in Selim's code
   * @see #GABOR_DATA_FEED
   */
  public final static int GABOR_USE_HSI_AVERAGE = 1;
  /**
   * Feed Gabor feature extraction processor with a mistaken data format:
   * luminance data is first extracteda and then that data will be treated as
   * blue component of a rgb color and fed into the average of H,S,I procedure.
   * @see #GABOR_DATA_FEED
   */
  public final static int GABOR_USE_MISTAKE = 2;

  /**
   * How the rgb colors will be transformed before fed into Gabor process.
   */
  public static int GABOR_DATA_FEED = GABOR_USE_HSI_AVERAGE;

  public static boolean DEMO_FILTERED = true;

  public static int IMAGE_RESIZED = 256;

  public static String RESIZED_FILE_SUFFIX = ".s";

  public static int BELL_RESAMPLER = 0;
  public static int BSPLINE_RESAMPLER = 1;
  public static int HERMITE_RESAMPLER = 2;
  public static int LANCZOS3_RESAMPLER = 3;
  public static int LINEAR_RESAMPLER = 4;
  public static int MITCHELL_RESAMPLER = 5;
  public static int NEAREST_RESAMPLER = 6;
  public static int SIMPLE_RESAMPLER = 7;

  public static int RESAMPLER = NEAREST_RESAMPLER;

  int[][] getColorHistogram(int[] rgbs, int[] labels, int label_num) {
    int[][] hists = new int[label_num][SmithColorQuantization.MAX_LABEL + 1];
    for (int i=0; i<rgbs.length; i++) {
      int bin = SmithColorQuantization.quantize(rgbs[i]);
      hists[labels[i]][bin] ++;
    }
    return hists;
  }

  static int[] getLuminance(int[] pixels, int width, int height) {
    return CYMaker.getLuminance(pixels, new Dimension(width, height));
  }

  static int[] getHSIAverage(int[] pixels, int width, int height) {
    int[][] hsi = CYMaker.getHSI(pixels, new Dimension(width, height));
    int[] average = new int[pixels.length];
    for (int i=0; i<average.length; i++) {
      average[i] = (int)Math.round((hsi[i][0] + hsi[i][1] + hsi[i][2]) / 3.0);
    }
    return average;
  }

  static int[] getMistaken(int[] pixels, int width, int height) {
    pixels = getLuminance(pixels, width, height);
    return getHSIAverage(pixels, width, height);
  }

  public static void main(String[] args) {
    // init configuration
    Configuration.init(args);

    // get the image name
    String img_name = ReadArgu.getOption('i', args);
    if (img_name == null) {
      System.err.println("Usage: -i image_name");
      System.exit(1);
    }

    // whether scale the image first
    String scale_str = ReadArgu.getOption('s', args);

    // filter the results with a high-D range
    String range_str = ReadArgu.getOption('r', args);

    // read in image
    Image image = null;
    Dimension d = new Dimension();
    int[] pixels = null;
    try {
      image = yi.util.io.ImageReader.readImageFile(img_name);
      pixels = yi.util.io.ImageReader.getPixels(image, d);
    }
    catch (java.io.IOException e) {
      System.err.println("Can't read image " + img_name);
      e.printStackTrace();
      System.exit(2);
    }

    // process
    try {

      if (scale_str != null) {
        int new_size = Integer.parseInt(scale_str);
        ResampledDescriptor resampled = Resampler.scaleDown(pixels, d, (new Gabor()).RESAMPLER, new_size);
        if (resampled != null) {
          pixels = resampled.m_resampled;
          d = resampled.m_newDimension;
        }
        img_name = img_name + Resampler.RESIZED_FILE_SUFFIX + new_size;
      }

      Configuration.setParameter("SOURCE_NAME", img_name);
      yi.util.ImageUtil.setPixels(pixels);
      yi.util.ImageUtil.setDimension(d);

      long startTime = System.currentTimeMillis();

      // transform the 3D rgb colors to acceptable 1D data to gabor process
      int[] pixels_1D;
      if (GABOR_DATA_FEED == GABOR_USE_LUMINANCE) {
        pixels_1D = getLuminance(pixels, d.width, d.height);
      }
      else if (GABOR_DATA_FEED == GABOR_USE_HSI_AVERAGE) {
        pixels_1D = getHSIAverage(pixels, d.width, d.height);
      }
      else {  // GABOR_DATA_FEED == GABOR_USE_MISTAKE
        pixels_1D =getMistaken(pixels, d.width, d.height);
      }

      // prepare the data to feed into Gabor: change it to matrix and change the
      // data type from integer to double
      double [][] data_2D = new double[d.height][d.width];
      int count = 0;
      for (int y=0; y<d.height; y++) {
        for (int x=0; x<d.width; x++, count++) {
          data_2D[y][x] = pixels_1D[count];
        }
      }
      DoubleMatrix to_filter = new DoubleMatrix(d.height, d.width, data_2D);

      // compute texture features

  //    // Use yi's code
  //    Gabor.numOrients = 4; Gabor.numScales = 4;
  //    DoubleMatrix features = Gabor.getFiltered(image, Gabor.numScales, Gabor.numOrients);

      // Use Selim's code
      DoubleMatrix features = (new Gabor()).getReusableGaborRst(to_filter, Gabor.numScales, Gabor.numOrients);

      long endTime = System.currentTimeMillis();
      System.out.println("running time: " + (endTime-startTime));

      // demo the results
      if (DEMO_FILTERED) {

        int[][] int_features_data = new int[ features.height ][ features.width ];

        for( int c = 0; c < features.width; c++ ) {
          for(int r = 0; r < features.height; r++ ) {
              int_features_data[r][c] = (int) Math.min(255, Math.round(features.data[r][c] * 50));
            }
        }

        for (int c=0; c<features.width; c++) {
          int[] data = new int[features.height];
          for (int r= 0; r < features.height; r++ ) {
            data[r] = int_features_data[r][c];
          }
          yi.util.ImageUtil.makeImageFromGray(data);
          OutputResult.outputGIFImage(img_name + ".g" + Conversion.toString(c, 2) + ".gif",data, d.width, d.height);
        }
      }

      // filter using range
      if (range_str != null) {

        // read in the ranges
        String[] range_strs = yi.util.io.ReadArgu.readStrings(range_str);
        double[][] ranges = new double[features.width][2];
        for (int c=0; c<features.width; c++) {
          StringTokenizer tokenizer = new StringTokenizer(range_strs[c], ", \t");
          ranges[c][0] = Double.parseDouble(tokenizer.nextToken());
          ranges[c][1] = Double.parseDouble(tokenizer.nextToken());
        }

        // filtering the image
        int[] filtered = new int[features.height];
        for(int r = 0; r < features.height; r++ ) {
          filtered[r] = 0;
          for( int c = 0; c < features.width; c++ ) {
            if (features.data[r][c] < ranges[c][0]) break;
            if (features.data[r][c] >= ranges[c][1]) break;
            filtered[r] = 1;
          }
        }

        // output the filtered image
        yi.util.ImageUtil.makeImageFromBinary(filtered);
        OutputResult.outputGIFImage(img_name + ".filtered.gif", filtered, d.width, d.height);
      }

      // output a portion of the image
      int x1 = 20, y1=105, x2=86, y2=139;
      double[][] portion = new double[(x2-x1)*(y2-y1)][];
      int count1 = 0, count2 = 0;
      for (int y=0; y<d.height; y++) {
        for (int x=0; x<d.width; x++, count2++) {
          if (x>=x1 && x<x2 && y>=y1 && y<y2) {
            portion[count1++] = features.data[count2];
          }
        }
      }
      System.out.println("count = " + count1);
      OutputResult.print("regions", portion, portion[0].length, portion.length);

    } catch (Exception e) {
      System.err.println("broken when processing " + img_name);
      e.printStackTrace();
    }

    System.exit(0);
  }

}
