package ij.process;
import ij.measure.*;
import java.awt.*;

/** Statistics, including the histogram, of an image or selection. */
public class ImageStatistics implements Measurements {

	/** Use the hIstogram() method to get the histogram as a double  array. */
	public int[] histogram;
	
	/** Int pixel count (limited to 2^31-1) */
	public int pixelCount;
	/** Long pixel count */
	public long longPixelCount;	
	
	/** Mode */
	public double dmode;
	/** Mode of 256 bin histogram (counts limited to 2^31-1) */	
	public int mode;
		
	public double area;
	public double min;
	public double max;
	public double mean;
	public double median;
	public double stdDev;
	public double skewness;
	public double kurtosis;
	public double xCentroid;
	public double yCentroid;
	public double xCenterOfMass;
	public double yCenterOfMass;
	public double roiX, roiY, roiWidth, roiHeight;
	/** Uncalibrated mean */
	public double umean;
	/** Length of major axis of fitted ellipse */
	public double major;
	/** Length of minor axis of fitted ellipse */
	public double minor;
	/** Angle in degrees of fitted ellipse */
	public double angle;
	/** Bin width 1 histogram of 16-bit images */
	public int[] histogram16;
	/** Long histogram; use getHIstogram() to retrieve. */
	protected long[] longHistogram;
	public double areaFraction;
	/** Used internally by AnalyzeParticles */
	public int xstart, ystart;
	/** Used by HistogramWindow */
	public boolean stackStatistics;
	/** Minimum threshold when "Limit to threshold" enabled */
	public double lowerThreshold = Double.NaN;	
	/** Maximum threshold when "Limit to threshold" enabled */
	public double upperThreshold = Double.NaN;	
	public double histMin;
	public double histMax;
	public int histYMax;
	public int maxCount;
	public int nBins = 256;
	public double binSize = 1.0;
	
	protected int width, height;
	protected int rx, ry, rw, rh;
	protected double pw, ph;
	protected Calibration cal;
	
	EllipseFitter ef;

	
	/** Calculates and returns uncalibrated (raw) statistics for the
	 * specified image, including histogram, area, mean, min and
	 * max, standard deviation and mode.  Use ImageProcessor.setRoi(x,y,width,height)
	 * to limit statistics to a rectangular area and ImageProcessor.setRoi(Roi)
	 * to limit to a non-rectangular area.
	 * @see ij.process.ImageProcessor#setRoi(int,int,int,int)
	 * @see ij.process.ImageProcessor#setRoi(Roi)
	 * @see ij.process.ImageProcessor#getStats
	*/
	public static ImageStatistics getStatistics(ImageProcessor ip) {
		return getStatistics(ip, AREA+MEAN+STD_DEV+MODE+MIN_MAX+RECT, null);
	}

	/** Calculates and returns statistics for the specified
	 * image using the specified measurent options
	 * and calibration. Use ImageProcessor.setRoi(x,y,width,height)
	 * to limit statistics to a rectangular area and ImageProcessor.setRoi(Roi)
	 * to limit to a non-rectangular area.
	 * @see ij.process.ImageProcessor#setRoi(int,int,int,int)
	 * @see ij.process.ImageProcessor#setRoi(Roi)
	 * @see ij.measure.Measurements
	*/
	public static ImageStatistics getStatistics(ImageProcessor ip, int mOptions, Calibration cal) {
		Object pixels = ip.getPixels();
		if (pixels instanceof byte[])
			return new ByteStatistics(ip, mOptions, cal);
		else if (pixels instanceof short[])
			return new ShortStatistics(ip, mOptions, cal);
		else if (pixels instanceof int[])
			return new ColorStatistics(ip, mOptions, cal);
		else if (pixels instanceof float[])
			return new FloatStatistics(ip, mOptions, cal);
		else
			throw new IllegalArgumentException("Pixels are not byte, short, int or float");
	}

	void getRawMinAndMax(int minThreshold, int maxThreshold) {
		int min = minThreshold;
		while ((histogram[min] == 0) && (min < 255))
			min++;
		this.min = min;
		int max = maxThreshold;
		while ((histogram[max] == 0) && (max > 0))
			max--;
		this.max = max;
	}

	void getRawStatistics(int minThreshold, int maxThreshold) {
		int count;
		double value;
		double sum = 0.0;
		double sum2 = 0.0;
		
		for (int i=minThreshold; i<=maxThreshold; i++) {
			count = histogram[i];
			longPixelCount += count;
			sum += (double)i*count;
			value = i;
			sum2 += (value*value)*count;
			if (count>maxCount) {
				maxCount = count;
				mode = i;
			}
		}
		pixelCount = (int)longPixelCount;
		area = longPixelCount*pw*ph;
		mean = sum/longPixelCount;
		umean = mean;
		dmode = mode;
		calculateStdDev(longPixelCount, sum, sum2);
		histMin = 0.0;
		histMax = 255.0;
	}
	
	void calculateStdDev(double n, double sum, double sum2) {
		if (n>0.0) {
			stdDev = (n*sum2-sum*sum)/n;
			if (stdDev>0.0)
				stdDev = Math.sqrt(stdDev/(n-1.0));
			else
				stdDev = 0.0;
		} else
			stdDev = 0.0;
	}
		
	void setup(ImageProcessor ip, Calibration cal) {
		width = ip.getWidth();
		height = ip.getHeight();
		this.cal = cal;
		Rectangle roi = ip.getRoi();
		if (roi != null) {
			rx = roi.x;
			ry = roi.y;
			rw = roi.width;
			rh = roi.height;
		}
		else {
			rx = 0;
			ry = 0;
			rw = width;
			rh = height;
		}
		
		if (cal!=null) {
			pw = cal.pixelWidth;
			ph = cal.pixelHeight;
		} else {
			pw = 1.0;
			ph = 1.0;
		}
		
		roiX = cal!=null?cal.getX(rx):rx;
		roiY = cal!=null?cal.getY(ry, height):ry;
		roiWidth = rw*pw;
		roiHeight = rh*ph;
	}
	
	void getCentroid(ImageProcessor ip) {
		byte[] mask = ip.getMaskArray();
		int count=0, mi;
		double xsum=0.0, ysum=0.0;
		for (int y=ry,my=0; y<(ry+rh); y++,my++) {
			mi = my*rw;
			for (int x=rx; x<(rx+rw); x++) {
				if (mask==null||mask[mi++]!=0) {
					count++;
					xsum += x;
					ysum += y;
				}
			}
		}
		xCentroid = xsum/count+0.5;
		yCentroid = ysum/count+0.5;
		if (cal!=null) {
			xCentroid = cal.getX(xCentroid);
			yCentroid = cal.getY(yCentroid, height);
		}
	}
	
	void fitEllipse(ImageProcessor ip, int mOptions) {
		ImageProcessor originalMask = null;
		boolean limitToThreshold = (mOptions&LIMIT)!=0 && ip.isThreshold();
		if (limitToThreshold) {
			ImageProcessor mask = ip.getMask();
			Rectangle r = ip.getRoi();
			if (mask==null) {
				mask = new ByteProcessor(r.width, r.height);
				mask.invert();
			} else {
				originalMask = mask;
				mask = mask.duplicate();
			}
			int n = r.width*r.height;
			double t1 = ip.getMinThreshold();
			double t2 = ip.getMaxThreshold();
			double value;
			for (int y=0; y<r.height; y++) {
				for (int x=0; x<r.width; x++) {
					value = ip.getf(r.x+x, r.y+y);
					if (value<t1 || value>t2)
						mask.setf(x, y, 0f);
				}
			}
			ip.setMask(mask);
		}
		if (ef==null)
			ef = new EllipseFitter();
		ef.fit(ip, this);
		if (limitToThreshold) {
			if (originalMask==null)
				ip.setMask(null);
			else
				ip.setMask(originalMask);
		}
		double psize = (Math.abs(pw-ph)/pw)<.01?pw:0.0;
		major = ef.major*psize;
		minor = ef.minor*psize;
		angle = ef.angle;
		xCentroid = ef.xCenter;
		yCentroid = ef.yCenter;
		if (cal!=null) {
			xCentroid = cal.getX(xCentroid);
			yCentroid = cal.getY(yCentroid, height);
		}
	}
	
	public void drawEllipse(ImageProcessor ip) {
		if (ef!=null)
			ef.drawEllipse(ip);
	}
	
	void calculateMedian(int[] hist, int first, int last, Calibration cal) {
		//ij.IJ.log("calculateMedian: "+first+"  "+last+"  "+hist.length+"  "+pixelCount);
		if (pixelCount==0 || first<0 || last>hist.length) {
			median = Double.NaN;
			return;
		}
		double sum = 0;
		int i = first-1;
		double halfCount = pixelCount/2.0;
		do {
			sum += hist[++i];
		} while (sum<=halfCount && i<last);
		median = cal!=null?cal.getCValue(i):i;
	}
	
	void calculateAreaFraction(ImageProcessor ip, int[] hist) {
		int sum = 0;
		int total = 0;
		int t1 = (int)Math.round(ip.getMinThreshold());
		int t2 = (int)Math.round(ip.getMaxThreshold());
		if (t1==ImageProcessor.NO_THRESHOLD) {
			for (int i=0; i<hist.length; i++)
				total += hist[i];
			sum = total - hist[0];
		} else {
			for (int i=0; i<hist.length; i++) {
				if (i>=t1 && i<=t2)
					sum += hist[i];
				total += hist[i];
			}
		}
		areaFraction = sum*100.0/total;
	}
	
	/** Returns the histogram as an array of doubles. */
	public double[] histogram() {
		double[] hist = new double[histogram.length];
		for (int i=0; i<hist.length; i++) {
			if (longHistogram!=null)
				hist[i] = longHistogram[i];
			else
				hist[i] = histogram[i];
		}
		return hist;
	}

	/** Returns the histogram as an array of longs. */
	public long[] getHistogram() {
		double[] hist = histogram();
		long[] hist2 = new long[hist.length];
		for (int i=0; i<hist.length; i++)
			hist2[i] = (long)hist[i];
		return hist2;
	}
	
	public String toString() {
		return "stats[count="+pixelCount+", mean="+mean+", min="+min+", max="+max+"]";
	}
	
	protected void saveThreshold(double minThreshold, double maxThreshold, Calibration cal) {
		if (cal!=null) {
			minThreshold = cal.getCValue(minThreshold);
			maxThreshold = cal.getCValue(maxThreshold);
		}
		lowerThreshold = minThreshold;
		upperThreshold = maxThreshold;
	}

}