package yi.iu.region;

import java.awt.*;
import java.awt.image.*;
import java.util.*;
import yi.util.io.*;
import yi.util.data.*;
import yi.graph.polygon.*;
import yi.iu.imgdb.learn.*;
import yi.iu.imgdb.learn.absrgn.*;

/**
 * Title:        Region
 * Description:  Contains stastistics calculation function of regions
 * Copyright:    Copyright (c) 2002
 * Company:      University of Washington
 * @author Yi Li
 * @version 1.0
 */

public class Region implements java.io.Serializable, AbstractRegion{
    static final float SMALLTOZERO = 0.01f;   // some float calcuation results no
                                        // zero, but a very small float

    public int    m_size;
    public float  m_centroidY;
    public float  m_centroidX;
    public double m_majAngle,m_majLength;  // major axis (total lenght, not the half!)
    public double m_minAngle,m_minLength;

    public float  m_momentY;     // second-order row moment
    public float  m_momentX;     // second-order column moment
    public float  m_momentXY;    // second-order mixed moment
    public float  m_d;           // R = {(r,c) | dr^2 + 2erc + fc^2 <= 1}
    public float  m_e;
    public float  m_f;

    /** Normalized region properties */
    public int m_imgWidth, m_imgHeight; // kept to normalize the parameters

    static final long serialVersionUID =  -4562635294352193291L;

    public Region() {
    }

    public String toString() {
      return new String(
        "<Region" +
          " size=\"" + m_size + "\"" +
          " centroidX=\"" + m_centroidX +  "\"" +
          " centroidY=\"" + m_centroidY +  "\"" +
          " momentX=\"" + m_momentX +  "\"" +
          " momentY=\"" + m_momentY +  "\"" +
          " momentXY=\"" + m_momentXY +  "\"" +
          " d=\"" + m_d +  "\"" +
          " e=\"" + m_e +  "\"" +
          " f=\"" + m_f +  "\"" +
          " majAngle=\"" + m_majAngle +  "\"" +
          " majLength=\"" + m_majLength +  "\"" +
          " minAngle=\"" + m_minAngle +  "\"" +
          " minLength=\"" + m_minLength +  "\"" +
          " imgWidth=\"" + m_imgWidth +  "\"" +
          " imgHeight=\"" + m_imgHeight +  "\"" +
        "/>\n");
    }

    public static Region combine(Region r1, Region r2) {
      Region r = new Region();
      r.m_centroidX = (r1.m_centroidX * r1.m_size + r2.m_centroidX * r2.m_size) / (r1.m_size + r2.m_size);
      r.m_centroidY = (r1.m_centroidY * r1.m_size + r2.m_centroidY * r2.m_size) / (r1.m_size + r2.m_size);
      double dx = r2.m_centroidX - r1.m_centroidX, dy = r2.m_centroidY -r1.m_centroidY;
      r.m_majAngle = Math.atan2(dy, dx);
      r.m_minAngle = rotateHalfPI(r.m_majAngle);
      double radius1 = Math.sqrt(r1.m_majLength*r1.m_minLength)/2, radius2 = Math.sqrt(r2.m_majLength*r2.m_minLength)/2;
      double distance = Math.sqrt(dx * dx + dy * dy);
      r.m_majLength = distance + radius1 + radius2;
      r.m_minLength = Math.max(radius1, radius2) * 2;
      if (r.m_majLength < r.m_minLength)  {  // exchange maj and min
        double temp = r.m_majLength;
        r.m_majLength = r.m_minLength; r.m_minLength = temp;
        temp = r.m_minAngle;
        r.m_minAngle = r.m_majAngle;
        r.m_majAngle = temp;
      }
      r.m_size = (int)Math.round(Math.PI * r.m_majLength * r.m_minLength / 4);

      if (r1.m_imgHeight!=r2.m_imgHeight || r1.m_imgWidth!=r2.m_imgWidth) {
        System.err.println("regions to combined are not based on the same size image");
      }
      r.m_imgWidth = r1.m_imgWidth;
      r.m_imgHeight = r1.m_imgHeight;

      return r;
    }

    public void scale(double scale) {
      m_size = (int)Math.round(m_size * scale * scale);
      m_centroidY = (float)(m_centroidY * scale);
      m_centroidX = (float)(m_centroidX * scale);
      m_majLength = m_majLength * scale;
      m_minLength = m_minLength * scale;

      /** @todo  m_momentY, m_momentX, m_momentXY, m_d, m_e, m_f are not scaled */

      m_imgWidth = (int)Math.round(m_imgWidth * scale);
      m_imgHeight = (int)Math.round(m_imgHeight * scale);
    }

    public void shift(double dx, double dy) {
      m_centroidY = (float)(m_centroidY + dy);
      m_centroidX = (float)(m_centroidX + dx);
    }

    public double[] normalize() {
      int size = m_imgWidth * m_imgHeight;
      double m_sizeNorm = m_size/(double)size;
      double m_centroidYNorm = (double)m_centroidY/m_imgHeight;
      double m_centroidXNorm = (double)m_centroidX/m_imgWidth;
      double m_majLengthNorm = normalizeLineLength(m_majAngle, m_majLength);
      double m_minLengthNorm = normalizeLineLength(m_minAngle, m_minLength);
      double m_momentXNorm = (double)m_momentX/m_imgWidth/m_imgWidth;
      double m_momentYNorm = (double)m_momentY/m_imgHeight/m_imgHeight;
      double m_momentXYNorm = (double)m_momentXY/m_imgWidth/m_imgHeight;
      return new double[] {
        m_sizeNorm,
        m_centroidXNorm,
        m_centroidYNorm,
        m_majLengthNorm,
        m_minLengthNorm,
        m_momentXNorm,
        m_momentYNorm,
        m_momentXYNorm
      };
    }

    public double getNormalizedX() { return (double)m_centroidX/m_imgWidth; }
    public double getNormalizedY() { return (double)m_centroidY/m_imgHeight; }
    public double getNormalizedSize() { return m_size/(double)(m_imgWidth * m_imgHeight); }

    public double normalizeLineLength(double angle, double length) {
      double dw = length * Math.cos(angle) / m_imgWidth;
      double dh = length * Math.sin(angle) / m_imgHeight;
      return Math.sqrt(dw*dw+dh*dh);
    }

    public static class TestCombine {
      public static void main(String[] args) {
        Region r1= new Region();
        r1.m_centroidX = 100;
        r1.m_centroidY = 100;
        r1.m_majAngle = 0;
        r1.m_majLength = 50;
        r1.m_minAngle = Math.PI/2;
        r1.m_minLength = 50;
        r1.m_size = (int)Math.round(Math.PI*r1.m_majLength*r1.m_minLength/4);


        Region r2= new Region();
        r2.m_centroidX = 300;
        r2.m_centroidY = 100;
        r2.m_majAngle = 0;
        r2.m_majLength = 100;
        r2.m_minAngle = Math.PI/2;
        r2.m_minLength = 50;
        r2.m_size = (int)Math.round(Math.PI*r2.m_majLength*r2.m_minLength/4);

        Region r = Region.combine(r1, r2);
        double[] norms = r.normalize();

      }
    }

    public static class TestIntersect {
      public static void main(String[] args) {
        Region r1= new Region();
        r1.m_centroidX = 100;
        r1.m_centroidY = 100;
        r1.m_majAngle = 0;
        r1.m_majLength = 50;
        r1.m_minAngle = Math.PI/2;
        r1.m_minLength = 50;
        r1.m_size = (int)Math.round(Math.PI*r1.m_majLength*r1.m_minLength/4);


        Region r2= new Region();
        r2.m_centroidX = 120;
        r2.m_centroidY = 100;
        r2.m_majAngle = 0;
        r2.m_majLength = 100;
        r2.m_minAngle = Math.PI/2;
        r2.m_minLength = 50;
        r2.m_size = (int)Math.round(Math.PI*r2.m_majLength*r2.m_minLength/4);

        Region r = (Region)r1.intersect(r2);
        int width = 200, height = 200;
        int[] pixels = new int[width * height];
        Arrays.fill(pixels, 0xffffffff);
        r1.drawRegion(pixels, width, height, 0xffff0000);
        r2.drawRegion(pixels, width, height, 0xff00ff00);
        r.drawRegion(pixels, width, height, 0xff0000ff);

        try {
          yi.util.io.ImageIOUtil.write("intersect.png", "png", pixels, width,
                                       height);
        } catch (java.io.IOException e) {
          e.printStackTrace();
        }
      }
    }

    /**
     * While passed in angel is (-PI,PI], the returned result is also in this range
     * @param angle
     * @return
     */
    static double rotateHalfPI(double angle) {
      double rst = angle + Math.PI/2;
      if (rst > Math.PI) rst -= Math.PI * 2;
      return rst;
    }

    /**
     * Given a label array, return the region informations
     * @param labels
     * @param maxLabel inclusive
     * @param d
     * @return
     */
    public static Region[] getRegions(int[] labels, int maxLabel, Dimension d) {
        Region[] r = new Region[maxLabel+1];

        // init the Regions
        for (int i=0; i<=maxLabel; i++)
            r[i] = new Region();

        // ==============
        // area size
        for(int count=0; count<labels.length; count++)
        {
            if(labels[count]>maxLabel)
            {
                // bad values passed in
                System.out.println("Region.getRegions did not complete");
                System.out.println("maxLabel parameter incorrect");
                break;
            }
            r[labels[count]].m_size ++;
        }

        // ==============
        // centroid
        int[] yArray = new int[maxLabel+1];
        int[] xArray = new int[maxLabel+1];

        // now collect sums of row numbers and sums of column numbers
        // for each region label.
        for(int count=0; count<labels.length; count++)
        {
            yArray[labels[count]] += count/d.width;
            xArray[labels[count]] += count%d.width;
        }

        // now compute the centroid for each label
        // note, although the centroid for 0 is computed, it is not a label
        for(int count=0; count<=maxLabel; count++)
        {
            r[count].m_centroidY=(float)yArray[count]/(float)r[count].m_size;
            r[count].m_centroidX=(float)xArray[count]/(float)r[count].m_size;
        }

        // ==========================
        // The second spatial moments

        // first calculate the sum of the squared dev. from the mean
        float[] sumSquaredR = new float[maxLabel+1];
        float[] sumSquaredC = new float[maxLabel+1];
        float[] sumSquaredM = new float[maxLabel+1];
        for(int count=0; count<labels.length; count++)
        {
            sumSquaredR[labels[count]]+=Math.pow(
                (count/d.width-r[labels[count]].m_centroidY),2);
            sumSquaredC[labels[count]]+=Math.pow(
                (count%d.width-r[labels[count]].m_centroidX),2);
            sumSquaredM[labels[count]]+=
                ((count/d.width-r[labels[count]].m_centroidY)*
                (count%d.width-r[labels[count]].m_centroidX));
        }

        // now find the average of the second order deviations
        for(int count=0; count<=maxLabel; count++)
        {
            r[count].m_momentY  = sumSquaredR[count]/r[count].m_size;
            r[count].m_momentX  = sumSquaredC[count]/r[count].m_size;
            r[count].m_momentXY = sumSquaredM[count]/r[count].m_size;
        }

        // ====================
        // elipses

        // This function determines the coefficients of the equation
        // R={(r,c)|dr^2+2erc+fc^2<=1} which determines the elipse that
        // best approximates a region.

        for(int count=0; count<=maxLabel; count++)
        {
            float matrixCoef = (float)(1/(4*(r[count].m_momentY*
                r[count].m_momentX-Math.pow(r[count].m_momentXY,2))));

            r[count].m_d = matrixCoef*r[count].m_momentX;
            r[count].m_e = matrixCoef*(-r[count].m_momentXY);
            r[count].m_f = matrixCoef*r[count].m_momentY;
        }

        // ==============
        // axis

        // NOTE : this function is incomplete and still buggy!!!!!!!!!
        for(int count=0; count<=maxLabel; count++)
        {
//            System.out.println("moment cc x = " + r[count].m_momentX);
//            System.out.println("moment rr y = " + r[count].m_momentY);
//            System.out.println("moment rc xy= " + r[count].m_momentXY);
            if (Math.abs(r[count].m_momentXY)<SMALLTOZERO &&  // r[count].m_momentXY ==0
                r[count].m_momentY>r[count].m_momentX)
            {
                r[count].m_majAngle = -Math.PI/2;
                r[count].m_majLength = 4*Math.pow(r[count].m_momentY,.5);
                r[count].m_minAngle = 0;
                r[count].m_minLength = 4*Math.pow(r[count].m_momentX,.5);
            }
            else if(Math.abs(r[count].m_momentXY)<SMALLTOZERO && // // r[count].m_momentXY ==0
                    r[count].m_momentY <= r[count].m_momentX)
            {
                r[count].m_majAngle = 0;
                r[count].m_majLength = 4*Math.pow(r[count].m_momentX,.5);
                r[count].m_minAngle = -Math.PI/2;
                r[count].m_minLength = 4*Math.pow(r[count].m_momentY,.5);
            }
            else if(r[count].m_momentY <= r[count].m_momentX)   // r[count].m_momentXY!=0
            {
                r[count].m_majAngle=Math.atan(
                    (r[count].m_momentY-r[count].m_momentX+Math.pow(Math.pow(
                    (r[count].m_momentY-r[count].m_momentX),2)+4*Math.pow(
                    r[count].m_momentXY,2),.5))/(2.0*r[count].m_momentXY));
                r[count].m_majLength=Math.pow(8*(r[count].m_momentY+
                    r[count].m_momentX+Math.pow(Math.pow(
                    (r[count].m_momentY-r[count].m_momentX),2)+4*Math.pow(
                    r[count].m_momentXY,2),.5)),.5);
                r[count].m_minAngle=r[count].m_majAngle+Math.PI/2;
                r[count].m_minLength=Math.pow(8*(r[count].m_momentY+
                    r[count].m_momentX-Math.pow(Math.pow(
                    (r[count].m_momentY-r[count].m_momentX),2)+4*Math.pow(
                    r[count].m_momentXY,2),.5)),.5);
            }
            else  // r[count].m_momentXY!=0 && r[count].m_momentY>r[count].m_momentX
            {
               r[count].m_majAngle=Math.atan(
                    (r[count].m_momentY-r[count].m_momentX+Math.pow(Math.pow(
                    (r[count].m_momentY-r[count].m_momentX),2)+4*Math.pow(
                    r[count].m_momentXY,2),.5))/(2.0*r[count].m_momentXY));
                r[count].m_majLength=Math.pow(8*(r[count].m_momentY+
                    r[count].m_momentX+Math.pow(Math.pow(
                    (r[count].m_momentY-r[count].m_momentX),2)+4*Math.pow(
                    r[count].m_momentXY,2),.5)),.5);
                r[count].m_minAngle=r[count].m_majAngle+Math.PI/2;
                r[count].m_minLength=Math.pow(8*(r[count].m_momentY+
                    r[count].m_momentX-Math.pow(Math.pow(
                    (r[count].m_momentY-r[count].m_momentX),2)+4*Math.pow(
                    r[count].m_momentXY,2),.5)),.5);
            }
        }

        // the normalized properties
        for (int count=0; count<=maxLabel; count++) {
          r[count].m_imgHeight = d.height;
          r[count].m_imgWidth = d.width;
        }

        return r;
    }

    /**
     * Return the intersection of two ellipse.  Use approximate polygon to
     * calculate the intersection of two ellipse (not exact solution!)
     * @param ar2
     * @return
     */
    public AbstractRegion intersect(AbstractRegion ar2) {
      if (!(ar2 instanceof Region)) {
        System.err.println("Can't do intersect with other type of AbstractRegion with ellipse");
        return null;
      }
      // each PI/2, use n points to approximate
      int n = 4;
      double[] x1 = new double[n*4], y1 = new double[n*4];
      double[] x2 = new double[n*4], y2 = new double[n*4];
      getApproximatePoints(n, x1, y1);
      ((Region)ar2).getApproximatePoints(n, x2, y2);
      GenericPolygon poly1 = new GenericPolygon(new GenericContour(x1, y1));
      GenericPolygon poly2 = new GenericPolygon(new GenericContour(x2, y2));
      GenericPolygon poly3 = GenericPolygon.gpcPolygonClip(GenericPolygon.GPC_INT, poly1, poly2);

      if (poly3.getArea() == 0) return null;
      int left = (int)Math.floor(poly3.getLeft());
      int right = (int)Math.ceil(poly3.getRight());
      int top = (int)Math.floor(poly3.getTop());
      int bottom = (int)Math.ceil(poly3.getBottom());

      poly3.shift(-left, -top);
      int width = right - left + 1;
      int height = bottom - top + 1;

      int[] labels = new int[width * height]; // backgournd = 0
      poly3.filledDraw(labels, width, height, 1);
//      try {
//        int[] demo = yi.util.ImageUtil.makeNewImageFromLabel(labels, true);
//        yi.util.io.ImageIOUtil.write("poly_intersect.png", "png", demo, width,
//                                     height);
//      }catch (Exception e) {
//
//      }
      Region[] regions = yi.iu.region.Region.getRegions(labels, 1, new Dimension(width, height));
      Region int_region = regions[1];
      int_region.shift(left, top);
      int_region.m_imgHeight = m_imgHeight; int_region.m_imgWidth = m_imgWidth;

//      {
//        width = m_imgWidth; height = m_imgHeight;
//        int[] pixels = new int[width * height];
//        Arrays.fill(pixels, 0xffffffff);
//        drawRegion(pixels, width, height, 0xffff0000);
//        ((Region)ar2).drawRegion(pixels, width, height, 0xff00ff00);
//        int_region.drawRegion(pixels, width, height, 0xff0000ff);
//
//        try {
//          yi.util.io.ImageIOUtil.write("intersect.png", "png", pixels, width,
//                                       height);
//        }
//        catch (java.io.IOException e) {
//          e.printStackTrace();
//        }
//      }
      return int_region;
    }

    /**
     * Return the size of the ellipse
     * @return
     */
    public double size() {
      return m_size;
    }

    /**
     * Return the approximate points on the boundary of the ellipse
     * @param n the number of points in PI/2 (so the total returned points
     *          will be 4*n
     * @param x x coordinates of the returned points
     * @param y y coordinates of the returned points
     */
    public void getApproximatePoints(int n, double[] x, double[] y) {
      if (x.length < n*4 || y.length < n*4) {
        System.err.println("the x and y length must be 4*n");
        System.err.println("the Approximate Points are not calculated");
        return;
      }

      // the sin and cos within PI/2
      SinCosArray sincos = SinCosCache.getSinCosArray(n);
      double[] sin = sincos.sin;
      double[] cos = sincos.cos;

      double a = m_majLength / 2, b = m_minLength / 2;
      // the points on the ellipse centered at (0,0) and no rotation
      for (int i=0; i<=n; i++) {
        x[i] = a * cos[i]; y[i] = b * sin[i];
        x[2 * n - i] = -x[i]; y[2 * n - i] =  y[i];
        x[2 * n + i] = -x[i]; y[2 * n + i] = -y[i];
        if (i != 0) {
          x[4 * n - i] = x[i];
          y[4 * n - i] = -y[i];
        }
      }

      // rotate the points
      double theta = m_majAngle;
      double sin_theta = Math.sin(theta), cos_theta = Math.cos(theta);
      for (int i=0; i<4*n; i++) {
        double x_cp = x[i];
        x[i] = x[i] * cos_theta - y[i] * sin_theta;
        y[i] = x_cp * sin_theta + y[i] * cos_theta;
      }

      // shift the points
      for (int i=0; i<4*n; i++) {
        x[i] += m_centroidX;
        y[i] += m_centroidY;
      }
    }

    /**
     * Draw the ellipse contour using approximate points
     * @param pixels
     * @param width
     * @param height
     * @param color
     * @param width
     * @param n 1/4 of the points on the boundary
     */
    public void drawApproximateContour(int[] pixels, int width, int height,
                                       int color, int line_width, int n) {
      double[] x = new double[n*4];
      double[] y = new double[n*4];
      getApproximatePoints(n, x, y);
      GenericContour contour = new GenericContour(x, y);
      contour.draw(pixels, width, height, color, line_width);
    }

    /**
     * Mark the given image with the given regions.  The regions will be marked
     * in the image with red ellipse and red major/minor axis
     */
    public static void markRegions(int[] img, Dimension dim, Region[] regions)
    {
        int MARK = 0xffff0000;      // mark color = red

        // add elipses
        int numRegions = regions.length;
        for(int count=0; count<dim.width*dim.height; count++)
        {
            int y = count/dim.width;
            int x = count%dim.width;
                    for(int j=0; j<numRegions; j++)
                    {
                Region r = regions[j];
                        double check = r.m_d * Math.pow((y-r.m_centroidY),2)
                       + 2* r.m_e* (y-r.m_centroidY)* (x-r.m_centroidX)
                       + r.m_f * Math.pow((x-r.m_centroidX),2);
                        if(check<=1.0 && check>+.92)
                        {
                            img[count] = MARK;
                        }
                    }
                }

        BufferedImage bimg = new BufferedImage(dim.width, dim.height,
                      BufferedImage.TYPE_INT_ARGB);
        bimg.setRGB(0, 0, dim.width, dim.height, img, 0, dim.width);
        Graphics2D g = (Graphics2D)(bimg.getGraphics());

        // add axis
        for (int i=0; i<numRegions; i++) {
            Region r = regions[i];
//            System.out.println("major angel  = " + r.m_majAngle);
//            System.out.println("major Length = " + r.m_majLength);
//            System.out.println("minor angel  = " + r.m_minAngle);
//            System.out.println("minor Length = " + r.m_minLength);
//            System.out.println("Centroid x =" + r.m_centroidX);
//            System.out.println("Centroid y =" + r.m_centroidY);
            double dx1 = Math.cos(r.m_majAngle) * r.m_majLength / 2;
            double dy1 = Math.sin(r.m_majAngle) * r.m_majLength / 2;
            int x1 = (int)(r.m_centroidX + dx1 + 0.5);
            int x2 = (int)(r.m_centroidX - dx1 + 0.5);
            int y1 = (int)(r.m_centroidY + dy1 + 0.5);
            int y2 = (int)(r.m_centroidY - dy1 + 0.5);
            g.setColor(Color.red);
            g.drawLine(x1, y1, x2, y2);

            double dx2 = Math.cos(r.m_minAngle) * r.m_minLength / 2;
            double dy2 = Math.sin(r.m_minAngle) * r.m_minLength / 2;
            x1 = (int)(r.m_centroidX + dx2 + 0.5);
            x2 = (int)(r.m_centroidX - dx2 + 0.5);
            y1 = (int)(r.m_centroidY + dy2 + 0.5);
            y2 = (int)(r.m_centroidY - dy2 + 0.5);
            g.setColor(Color.blue);
            g.drawLine(x1, y1, x2, y2);
        }
        bimg.getRGB(0, 0, dim.width, dim.height, img, 0, dim.width);
//
//		// add axis
//		for(int count=0; count<dim.width*dim.height; count++)
//        {
//            int y = count/dim.width;
//            int x = count%dim.width;
//		    for(int j=0; j<numRegions; j++)
//		    {
//                Region r = regions[j];
//
//		        // check major axis first
//		        double slope;
//		        if(r.m_majAngle != Math.PI/2 && r.m_majAngle!=(-Math.PI/2))
//		        {
//		            slope = -Math.tan(r.m_majAngle);
//		        }
//		        else slope = 999;
//		        float deltaRow = y - r.m_centroidY;
//		        if(deltaRow == 0)
//		            deltaRow = (float).0001;
//		        float deltaColumn = x - r.m_centroidX;
//		        double distance = Math.sqrt(Math.pow(deltaColumn,
//		            2)+Math.pow(deltaRow,2));
//                if((slope>=0 &&distance<=.5*r.m_majLength&&
//		            slope>=.7*(deltaColumn/deltaRow)&&
//		            slope<=1.3*(deltaColumn/deltaRow)) ||
//		            (slope<0 &&distance<=.5*r.m_majLength&&
//		            slope<=.7*(deltaColumn/deltaRow)&&
//		            slope>=1.3*(deltaColumn/deltaRow)))
//		        {
//		            img[count]=MARK;
//		        }
//
//		        // now check minor axis
//		        if(r.m_minAngle!=Math.PI/2&&r.m_minAngle!=(-Math.PI/2))
//		        {
//		            slope = -Math.tan(r.m_minAngle);
//		        }
//		        //System.out.println(slope+"\n");
//		        if((slope>=0 &&distance<=.5*r.m_minLength&&
//		            slope>=.7*(deltaColumn/deltaRow)&&
//		            slope<=1.3*(deltaColumn/deltaRow)) ||
//		            (slope<0 &&distance<=.5*r.m_minLength&&
//		            slope<=.7*(deltaColumn/deltaRow)&&
//		            slope>=1.3*(deltaColumn/deltaRow)))
//		        {
//		            img[count]=MARK;
//		        }
//		    }
//		}
    }

    /**
     * Draw the region in the given image array
     * @param img the image array to draw.  If null, then an array will be
     *            constructed.
     * @param dim
     * @param r
     * @param Color
     * @param draw_axis draw the major and minor axis
     * @return
     */
    public static int[] drawRegion(int[] img, Dimension dim, Region r, int Color, boolean draw_axis) {
      if (img == null) {
        img = new int[dim.width*dim.height];
        java.util.Arrays.fill(img, 0xff000000);
      }

      for(int count=0; count<dim.width*dim.height; count++)
      {
          int y = count/dim.width;
          int x = count % dim.width;
          double check = r.m_d * Math.pow( (y - r.m_centroidY), 2)
              + 2 * r.m_e * (y - r.m_centroidY) * (x - r.m_centroidX)
              + r.m_f * Math.pow( (x - r.m_centroidX), 2);
          if (check <= 1.0) {
            img[count] = Color;
          }
      }
      // add axis
      if (draw_axis) {
        BufferedImage bimg = new BufferedImage(dim.width, dim.height,
                      BufferedImage.TYPE_INT_ARGB);
        bimg.setRGB(0, 0, dim.width, dim.height, img, 0, dim.width);
        Graphics2D g = (Graphics2D)(bimg.getGraphics());
        double dx1 = Math.cos(r.m_majAngle) * r.m_majLength / 2;
        double dy1 = Math.sin(r.m_majAngle) * r.m_majLength / 2;
        int x1 = (int)(r.m_centroidX + dx1 + 0.5);
        int y1 = (int)(r.m_centroidY + dy1 + 0.5);
//          int x2 = (int)(r.m_centroidX - dx1 + 0.5);
//          int y2 = (int)(r.m_centroidY - dy1 + 0.5);
        g.setColor(java.awt.Color.red);
        g.drawLine((int)r.m_centroidX, (int)r.m_centroidY, x1, y1);

        double dx2 = Math.cos(r.m_minAngle) * r.m_minLength / 2;
        double dy2 = Math.sin(r.m_minAngle) * r.m_minLength / 2;
        x1 = (int)(r.m_centroidX + dx2 + 0.5);
        y1 = (int)(r.m_centroidY + dy2 + 0.5);
//          x2 = (int)(r.m_centroidX - dx2 + 0.5);
//          y2 = (int)(r.m_centroidY - dy2 + 0.5);
        g.setColor(java.awt.Color.blue);
        g.drawLine((int)r.m_centroidX, (int)r.m_centroidY, x1, y1);
        bimg.getRGB(0, 0, dim.width, dim.height, img, 0, dim.width);
      }
      return img;
    }

    public static int[] drawRegions(int[] img, Dimension dim, Region r, int Color)
    {
      return drawRegion(img, dim, r, Color, true);
    }

    /**
     * Mark the given image with the given regions
     * @param img
     * @param dim
     * @param regions
     * @param colors
     * @param draw_axis
     * @return
     */
    public static int[] drawRegions(int[] img, Dimension dim, Region[] regions, int[] colors, boolean draw_axis)
    {
      if (img == null) {
          img = new int[dim.width*dim.height];
          java.util.Arrays.fill(img, 0xff000000);
      }

      for(int j=0; j<regions.length; j++) {
        drawRegion(img, dim, regions[j], colors[j], draw_axis);
      }

      return img;
    }

    public static int[] drawRegions(int[] img, Dimension dim, Region[] regions, int[] colors)
    {
      return drawRegions(img, dim, regions, colors, true);
    }

    public void drawRegion(int[] pixels, int width, int height, int color) {
      drawApproximateContour(pixels, width, height, color, 1, 8);

//      Region.drawRegion(pixels, new Dimension(width, height), this, color, false);
    }


    public static double[][] getMeans(int[] labels, double[][] features, int feature_num) {
      int max = yi.stat.StatUtil.maximum(labels);
//      System.out.println("region.max=" + max);
      double[][] means = new double[max+1][feature_num];
      int[] count = new int[max+1];

      // sum
      for (int i=0; i<labels.length; i++) {
        int index = labels[i];
        count[index]++;
        for (int j=0; j<feature_num; j++) {
          means[index][j] += features[i][j];
        }
      }

      // means
      for (int i=0; i<=max; i++) {
        for (int j=0; j<feature_num; j++) {
          means[i][j] /= count[i];
        }
      }
      return means;

    }

    public static double[][] getMedians(int[] labels, double[][] features, int feature_num) {
      int max = yi.stat.StatUtil.maximum(labels);
//      System.out.println("region.max=" + max);
      ArrayList[][] data = new ArrayList[max+1][feature_num];
      for (int i=0; i<= max; i++) {
        for (int j=0; j<feature_num; j++) {
          data[i][j] = new ArrayList();
        }
      }

      // put data to corresponding collection
      for (int i=0; i<labels.length; i++) {
        int index = labels[i];
        for (int j=0; j<feature_num; j++) {
          data[index][j].add(new Double(features[i][j]));
        }
      }

      // median
      double[][] medians = new double[max+1][feature_num];
      for (int i=0; i<=max; i++) {
        for (int j=0; j<feature_num; j++) {
          Collections.sort(data[i][j]);
          medians[i][j] = ((Double)data[i][j].get(data[i][j].size()/2)).doubleValue();
//          double[] region_data = Conversion.toDoubleArray(data[i][j]);
//          OutputResult.print("region_data." + Conversion.toString(i, 2) + "." +
//            Conversion.toString(j, 2), region_data);
//          System.out.println("region_data." + Conversion.toString(i, 2) + "." +
//            Conversion.toString(j, 2) + ": " + medians[i][j]);
        }
      }

      return medians;

    }

    /**
     * Output demonstration images with each one only shows one region.  The
     * image names will begin with the specified root_name and are followed by
     * ".region.#". # is the label number.
     * @param root_name
     * @param pixels
     * @param labels
     * @param label_num
     * @param width
     * @param heigth
     * @param background
     */
    public static void outputEachRegion(
      String root_name, int[] pixels, int[] labels, int label_num,
      int width, int height, int background) {
      int pnum = pixels.length;
      int[] demo_img = new int[pnum];

      for (int i=0; i<label_num; i++) {
        for (int j=0; j<pnum; j++) {
          demo_img[j] = (labels[j] == i) ? pixels[j] : background;
        }
        yi.util.io.OutputResult.outputColorImage(root_name + ".region." + i, demo_img, width, height, true);
      }
    }

    /**
     * Construct several ellipse regions and test whether the region methods
     * are right
     */
    public static class TestRegion {

      public static double distance(double x1, double y1, double x2, double y2) {
        double dx = x1 - x2, dy = y1 - y2;
        return Math.sqrt(dx * dx + dy * dy);
      }
      /**
       *
       */
      public static int ELLIPSE_NUM = 12;
      public static int c = 160;
      public static int a = 200;

      public static void main(String[] args) {
        // original point
        int original_x = (int)Math.round(a);
        int original_y = (int)Math.round(a);

        // image
        int image_w = original_x * 2;
        int image_h = original_y * 2;
        int[] labels = new int[image_w * image_h];

        for ( int i = 0; i < ELLIPSE_NUM; i ++ ) {
          double major_axis_angle = Math.PI * 2 / ELLIPSE_NUM * i;
          double c1_x = c * Math.cos(major_axis_angle);
          double c1_y = c * Math.sin(major_axis_angle);
          double c2_x = - c1_x, c2_y = - c1_y;

          // clear the labels
          java.util.Arrays.fill(labels, 0);

          // draw the ellipse in labels
          int count = 0;
          for (int y = 0; y < image_h; y++) {
            for (int x = 0; x < image_w; x++, count++) {
              int x_in_coord = x - original_x;
              int y_in_coord = y - original_y;
              if (distance(x_in_coord, y_in_coord, c1_x, c1_y) +
                  distance(x_in_coord, y_in_coord, c2_x, c2_y) <= 2 * a) {
                labels[count] = 1;
              }
            }
          }

          Region[] regions = Region.getRegions(labels, 1, new Dimension(image_w, image_h));
          Region ellipse_region = regions[1];

          // draw the labels
          int[] ellipse_img = yi.util.ImageUtil.makeNewImageFromLabel(labels);
          yi.util.io.OutputResult.outputGIFImage("ellipse." + i + ".gif", ellipse_img, image_w, image_h);

          // draw the region
          int[] region_img = new int[image_w * image_h];
          Region.drawRegion(region_img, new Dimension(image_w, image_h), ellipse_region, 0xffffffff, true);
          yi.util.io.OutputResult.outputGIFImage("region." + i + ".gif", region_img, image_w, image_h);
        }
        System.exit(0);

      }
    }


}



/**
 * Save the sin and cos of n angles (1/n of PI/2).  The total returned number is
 * (n+1): from 0, ..., PI/2
 */
class SinCosArray {
  double[] sin;
  double[] cos;

  /**
   * partition PI/2 to n segments and save the sin and cos results.  From 0,
   * , ..., (n-1)/n * (PI/2), PI/2.  n+1 in total
   * @param n
   */
  SinCosArray(int n) {
    sin = new double[n+1];
    cos = new double[n+1];
    double angle_n = Math.PI / 2 / n;
    for (int i=0; i<=n; i++) {
      double angle = i * angle_n;
      sin[i] = Math.sin(angle);
      cos[i] = Math.cos(angle);
    }
  }
}

/**
 * To save the results of "SinCosArray" so that they will be only calculated
 * once
 */
class SinCosCache {
  static HashMap m_sinCosArrays = new HashMap();

  /**
   * get the SinCosArray with parameter of n
   * @param n
   */
  public static SinCosArray getSinCosArray(int n) {
    Integer index = new Integer(n);
    if (!m_sinCosArrays.containsKey(index)) {
      m_sinCosArrays.put(index, new SinCosArray(n));
    }
    return (SinCosArray)m_sinCosArrays.get(index);
  }
}