import java.awt.*;
import java.io.*;
import ij.*;
import ij.gui.*;
import ij.process.*;
import ij.text.*;
import ij.plugin.PlugIn;
import ij.measure.*;
import ij.util.*;
import java.util.*;
import java.awt.event.*;


/*	This plugin correlates two 32-bit images or stacks.
	The plugin generates a scatter plot composed of the pixel values of Image1 (array1) and Image2 (array2)
	The statistics given are:
	- The correlation between Image1 and Image2 (R)
	- The slope of the regression line
	- The constant C intercepting the y-axis
	- The local area size can be predefined. An average value is calculated for each local area.
	- A 100 x100 pixels image yields an array of 10000 scatter plot values when defining a 1-pixel local area
		and 100 scatter plot values when defining a 10-pixels local area.
	- The max local area size is half the shortest image dimension (width or height)

	Gary Chinga 050102
	
	051206	Included a next option. When checked, all the slides in Stack1 will be correlated to the 
		current image of stack2. When correlating several slides in Stack1 to 1 image in Stack2 check the fix Slide option.
	060612	Included a next slide option. When checked, a given slide will be correlated with the next one in the same stack
		The fix slide option will then be inactivated.
	070408	Several changes to the dialog box
		Include the correlation of several statistics
		Option to create a correlation map
		Option to create the correlation arrays
	071218	Changes to the dialog box. 
		Implemented the increasing interval correlations
	080210	Fixed a small bud in the verification of stacks.	

*/

public class Image_CorrelationJ_1o implements PlugIn, ItemListener {

	private static int index1;
	private static int index2;
    	static int[] kernelX = new int[]{-1,0,1,-2,0,2,-1,0,1};
	static int[] kernelY = new int[]{1,2,1,0,0,0,-1,-2,-1};
	//private static boolean dR2=true, dEquation=true, dResults=true, fixSlide=false, nextSlide=false, cMap=true, cArray=false, cPlot=true;
	private static boolean  cMap=true, cArray=false, cPlot=true;
	private static int targetSource=0,fixSlide=1, nextSlide=2, iSequence=3;
	private ImageStack img1, img2;
	private int nSlices,dec,mSize,currentSlide=0;
	public static final int LINE=2,BOX=3, CIRCLE=0, CROSS=5,DOT=6,TRIANGLE=4, X=1;
	public static final String[] markers ={"X","Cross","Circle","Triangle","Box"};
	public static  int marker;

	public static final String[] cStats ={"Average","std","Skewness","kurtosis","Directionality"};
	public static  int cStat;

	public static final String[] cSequence ={"Target-Source","Fix slide (Source)","Next slide","Increasing sequence"};
	public static  int cSeq;

	
 	TextWindow tw,ctw;
	String title,headings,aLine ;
	String ctitle,cheadings,cLine ;
	String xAxis = "0,10";
  	String yAxis = "0,10";
	boolean xyFix=false;

	Checkbox cbFix = new Checkbox("Fix x- and y-axes",xyFix);
	Label lbxFix = new Label("x-axis (min,max): ");
	TextField tfxFix = new TextField("0,10",0);
	Label lbyFix = new Label("y-axis (min,max): ");
	TextField tfyFix = new TextField("0,10",0);
	Label lbNone1 = new Label("");
	Label lbNone2 = new Label("");
	Label lbNone3 = new Label("");
	Label lbNone4 = new Label("");
	Label lbNone5 = new Label("");
	Label lbNone7 = new Label("");
	Label Display = new Label("Display: ");
    public void run(String arg) {
        if (showDialog())
            plotCorrelation(img1, img2, currentSlide);
           //writeResults
    }

    public boolean showDialog() {
        int[] wList = WindowManager.getIDList();
        if (wList==null) {
            IJ.noImage();
            return false;
        }
        String[] titles = new String[wList.length];
        for (int i=0; i<wList.length; i++) {
            ImagePlus imp = WindowManager.getImage(wList[i]);
            if (imp!=null)
                titles[i] = imp.getTitle();
            else
                titles[i] = "";
        }
        if (index1>=titles.length)index1 = 0;
        if (index2>=titles.length)index2 = 0;
        GenericDialog gd = new GenericDialog("Image CorrelationJ 1o");

        gd.setLayout(new GridLayout(14,2,0,0));

        gd.addChoice("Target: ", titles, titles[index1]);
        gd.addChoice("Source: ", titles, titles[index2]);
        gd.addChoice("Correlation sequence: ", cSequence, cSequence[0]);
        
	//gd.addCheckbox("Fix slide (Target image)", fixSlide);
	//gd.addCheckbox("Correlate Next Slice", nextSlide);
	gd.add(lbNone2);	gd.add(lbNone7);
	gd.add(Display);//
	gd.addCheckbox("Correlation map", cMap);
        gd.addCheckbox("Correlation arrays", cArray);
	gd.addCheckbox("Correlation plot", cPlot);
	gd.addChoice("Statistic to correlate: ", cStats, cStats[0]);
	//gd.addNumericField("Local region size: ",3, 0);        
        gd.addNumericField("Local region size: ",3, 0);
        gd.addNumericField("Decimal places: ",2, 0);
	gd.add(cbFix);gd.add(lbNone1);
	gd.add(lbxFix); gd.add(tfxFix);
	gd.add(lbyFix); gd.add(tfyFix);
 	gd.addChoice("Markers: ", markers, markers[2]);
	lbxFix.setVisible(false);
	tfxFix.setVisible(false);
	lbyFix.setVisible(false);
	tfyFix.setVisible(false);

	cbFix.addItemListener(this);
        gd.add(lbNone3);//gd.add(lbNone4);

     	gd.showDialog();
        if (gd.wasCanceled())
            return false;
        index1 = gd.getNextChoiceIndex();
        index2 = gd.getNextChoiceIndex();
	cSeq = gd.getNextChoiceIndex();	
	cStat = gd.getNextChoiceIndex();
        marker = gd.getNextChoiceIndex();
	//fixSlide = gd.getNextBoolean();
	//nextSlide = gd.getNextBoolean();
	cMap=gd.getNextBoolean();
        cArray=gd.getNextBoolean();
        cPlot=gd.getNextBoolean();
        
	mSize = (int)gd.getNextNumber();
        dec = (int)gd.getNextNumber();

        xyFix = cbFix.getState();
        xAxis = tfxFix.getText();
        yAxis= tfyFix.getText();

        String title1 = titles[index1];
        String title2 = titles[index2];
        ImagePlus sliceimg1 = WindowManager.getImage(wList[index1]);
        ImagePlus sliceimg2 = WindowManager.getImage(wList[index2]);
	
	//if (nextSlide) fixSlide=false;	  
	
	if (cSeq==fixSlide) currentSlide = sliceimg1.getCurrentSlice();
        img1 = sliceimg1.getStack();
        img2 = sliceimg2.getStack();
      	
	if ((img1.getSize()!=img2.getSize())&&(cSeq!=fixSlide)) {
            IJ.showMessage("Both stacks must have the same number of slices or fix slide on source");
            return false;
        }
		int w1=0;

	if (img1.getWidth()!=img2.getWidth()||img1.getHeight()!=img2.getHeight()){
            IJ.showMessage("The stacks dimensions are not the same");
            return false;
        }
        else
		if (img1.getWidth()>img1.getHeight()) w1=img1.getHeight();
		else w1=img1.getWidth();

       if ((mSize<1)||(mSize>(w1/2))) {
            IJ.showMessage("The mask size must be larger than 0 and smaller than "+IJ.d2s((w1/2),0));
            return false;
        }

        return true;
   }

    public void plotCorrelation(ImageStack img1, ImageStack img2, int cSlide) {
        int nSlicesImg2 = img2.getSize();
	int nSlicesImg1 = img1.getSize();
	
     	int width = img1.getWidth();
        int height = img1.getHeight();
        ImageStack iStack = null;
	ImageStack cStack=null;
	ImageStack sStack=null;
	ImageStack tStack=null;
	
	ImageProcessor ipMap;
	ImageProcessor ipSource;
	ImageProcessor ipTarget;	
	ImageProcessor ipPlot;
	
	ImageProcessor ip1, ip2, plot=null;
	int xCells = (int)(width/mSize);
	int yCells = (int)(height/mSize);

	
	if((nSlicesImg1==1) && (cSeq==nextSlide)){
		IJ.showMessage("Stack required for correlating next slide...");
		return;
	}
	
	if (cSeq!=iSequence){
		if (cSeq==nextSlide) nSlicesImg1=nSlicesImg1-1;	
		for(int i=1; i<=nSlicesImg1; i++) {
			IJ.showProgress((double)(i)/nSlices);
			IJ.showStatus("a: "+(i)+"/"+nSlices);
			String label, label1, label2; //= img1.getSliceLabel(i)+","+img2.getSliceLabel(i);
				
			if	(img1.getSliceLabel(i)==null) label1="Target-"+IJ.d2s(i,0); else label1=img1.getSliceLabel(i);
			
			if (cSlide>0){	 //cSlide== currentSlice;cSlide!=0 means that a slide has been fixed.
				if	(img2.getSliceLabel(cSlide)==null) label2="Source-"+IJ.d2s(cSlide,0); else label2=img2.getSliceLabel(cSlide);
				ip2 = img2.getProcessor(cSlide); 
			}
			else{
				if (cSeq!=nextSlide){			//verifies if next slide correlation will be performed
					if	(img2.getSliceLabel(i)==null) label2="Source-"+IJ.d2s(i,0); else label2=img2.getSliceLabel(i);
					ip2 = img2.getProcessor(i); 
				}
				else{
					if	(img2.getSliceLabel(i+1)==null) label2="Source-"+IJ.d2s(i+1,0); else label2=img2.getSliceLabel(i+1);
					ip2 = img2.getProcessor(i+1); 
				}	
			}		
	
			label =label1+","+label2;	
			ip2 = ip2.convertToFloat();
				  
			ip1 = img1.getProcessor(i);     
			ip1 = ip1.convertToFloat();
			
			float[] ip1Array = getNewArray(ip1, mSize,cStat);
			float[] ip2Array = getNewArray(ip2, mSize,cStat); 
			double[] cStats= getCorrValues(ip1Array, ip2Array);
			writeResults(cStats,label1,label2);
			
			if (cPlot){
				ipPlot=generatePlot(ip1Array, ip2Array, label1, label2,cStats);
				if (iStack==null)
					iStack = new ImageStack(ipPlot.getWidth(), ipPlot.getHeight());
			
				iStack.addSlice(label, ipPlot);
			}		
			
			if (cMap){
				ipMap=generateMap(ip1Array, ip2Array,width,height, mSize, cStats, 2);
				if (cStack==null)
					cStack = new ImageStack(width, height);
				 cStack.addSlice(label,ipMap);
			}
			
			if (cArray){
				ipSource=generateMap(ip1Array, ip2Array,width,height, mSize, cStats, 1);
				if (sStack==null)
					sStack = new ImageStack(width, height);
				sStack.addSlice(label,ipSource);
	
				ipTarget=generateMap(ip1Array, ip2Array,width,height, mSize, cStats, 0);
				if (tStack==null)
					tStack = new ImageStack(width, height);
				tStack.addSlice(label,ipTarget);
			}
	
		}	
		if (cPlot){
			ImagePlus impCorrelation = new ImagePlus("Correlation plots",iStack);
			impCorrelation.show(); 
		}
		if (cMap) {
			ImagePlus impCorrelationMap = new ImagePlus("Correlation Maps",cStack);
			impCorrelationMap.show();IJ.selectWindow("Correlation Maps");IJ.run("Fire");
		}	
		if (cArray) {
			ImagePlus impTargetMap = new ImagePlus("Target Array",tStack);
			impTargetMap.show();IJ.selectWindow("Target Array");IJ.run("Fire");
			ImagePlus impSourceMap = new ImagePlus("Source Array",sStack);
			impSourceMap.show();IJ.selectWindow("Target Array");IJ.run("Fire");
		}		
	}
	else {
		
		if (nSlicesImg1>2){
			/*Frame interval 1 (1 sec): 1 vs 2, 2 vs 3, 3 vs 4, 4 vs 5
			Frame interval 2 (2 sec): 1 vs 3, 2 vs 4, 3 vs 5
			Frame interval 3 (3 sec): 1 vs 4, 2 vs 5
			Frame interval 4 (4 sec): 1 vs 5
			combine the R value from each interval size along with the 
			(sd of the R values for the interval and the size of the group (n)) and plot frame interval against R.
			*/
			boolean neg=false;
			float[] meanCorrelations = new float[nSlicesImg1-1];
			float[] stdCorrelations = new float[nSlicesImg1-1];
			float[] intervals = new float[nSlicesImg1-1];
			for (int ii=1; ii<=nSlicesImg1-1;ii++){
				float mCorrelation=0;//int cCounter=0;
				float[] mCorrArray = new float[nSlicesImg1-ii];
				for(int i=1; i<=nSlicesImg1-ii; i++) {
					IJ.showProgress((double)(i)/nSlices);
					IJ.showStatus("a: "+(i)+"/"+nSlices);					  
					ip1 = img1.getProcessor(i);  ip1 = ip1.convertToFloat();						
					ip2 = img2.getProcessor(i+ii); 	ip2 = ip2.convertToFloat();
					float[] ip1Array = getNewArray(ip1, mSize,cStat);
					float[] ip2Array = getNewArray(ip2, mSize,cStat); 
					double[] cStats= getCorrValues(ip1Array, ip2Array);
					mCorrArray[i-1]=(float)(cStats[2]);  //gets the R-value
					if (mCorrArray[i-1]<0) neg=true;
					//cCounter++;
				}
				intervals[ii-1]=ii;
				meanCorrelations[ii-1]=getMean(mCorrArray);
				stdCorrelations[ii-1]=getStd(meanCorrelations[ii-1],mCorrArray);
				writeResults(intervals[ii-1], meanCorrelations[ii-1],stdCorrelations[ii-1]);
				//IJ.log(IJ.d2s(ii,0)+"  "+IJ.d2s(meanCorrelations[ii-1],2));
			}	
			ipPlot=generatePlot(intervals, meanCorrelations,stdCorrelations,nSlicesImg1-1,neg);
			ImagePlus impPlot = new ImagePlus("Correlation curve",ipPlot);	
			impPlot.show();	
		}
		else {
			 IJ.showMessage("A stack with at least three slices is required");
			 return;
		}
	}
    }

	double[] getCorrValues(float[] a1, float[] a2){
		double[] d=new double[3];
		d[2]= getR(a1,a2);
		String sLabel="", R2label="";
		CurveFitter cf = new CurveFitter(Tools.toDouble(a1), Tools.toDouble(a2));
		cf.doFit(CurveFitter.STRAIGHT_LINE);
		double[] s = cf.getParams();
		d[0]=s[0];d[1]=s[1];
		return d;
	}


	ImageProcessor generateMap(float[] a1, float[] a2,int ww, int hh, int mSize, double[] corrStats, int map){
		int nY=(int)(hh/mSize);
		int nX=(int)(ww/mSize);

		float[] array= new float[nX*nY];
		double[] regMap=new double[ww*hh];
		for (int i=0;i<(ww*hh);i++)
			regMap[i]=0.0;
		FloatProcessor floatMap= new FloatProcessor(ww,hh,regMap);
		int counter=0, ind=0;
		for (int i = 0;i<nY;i++){
			int yStep=i*mSize;                                           
			double deltaRegValue=0;
			for (int ii = 0;ii<nX;ii++){
				int xStep=ii*mSize;
				float a1Value = a1[ind];
				float a2Value = a2[ind];
				if (map==0) deltaRegValue=a1Value;				
				if (map==1) deltaRegValue=a2Value;
				if (map==2) deltaRegValue=a2Value-((corrStats[0]+corrStats[1]*a1Value)); //Estimates the modelled reg value and the difference with y
				
				for (int y = yStep; y < (yStep+mSize); y++){
					for (int x = xStep; x < (xStep+mSize); x++){
						floatMap.putPixelValue(x, y, deltaRegValue);
					}
				}
				ind++;
			}
		}
		return floatMap;
	}

	ImageProcessor generatePlot(float[] cInt,float[] cValues,float[] sValues, int nSlices, boolean negR){
		//float[] a1=float[0][]
		Plot p = new Plot("Correlation plot", "Increasing intervals","Correlation (R)", cInt, cValues);
		if (negR) p.setLimits(1,nSlices,-1,1);
		else p.setLimits(1,nSlices,0,1);
		p.addPoints(cInt,cValues,PlotWindow.CIRCLE);
		p.addErrorBars(sValues);
		return p.getProcessor();
	}	
	
	ImageProcessor generatePlot(float[] a1, float[] a2, String lab1, String lab2, double[] corrStats){

		float R= (float)corrStats[2];
		String sLabel="", R2label="";
		double[] s =new double[2];
		s[0]=corrStats[0];
		s[1]=corrStats[1];		
		float[] px = new float[100];
		float[] py = new float[100];
		double[] a;
		double xmin,xmax,ymin,ymax;
		if (!xyFix){
			a = Tools.getMinMax(a1);
			xmin=a[0]; xmax=a[1];
			a = Tools.getMinMax(a2);
			ymin=a[0]; ymax=a[1];
		}
		else{
			int[] xLimit=s2ints(xAxis);
			int[] yLimit=s2ints(yAxis);
			xmin=xLimit[0];		xmax=xLimit[1];		ymin=yLimit[0];	ymax=yLimit[1];
		}
		double inc = (xmax-xmin)/99.0;
		double tmp = xmin;
		for (int ii=0; ii<100; ii++) {
			px[ii]=(float)tmp;
			tmp += inc;
		}
		for (int ii=0; ii<100; ii++)
			py[ii] = (float)CurveFitter.f(CurveFitter.STRAIGHT_LINE, s, px[ii]);
		a = Tools.getMinMax(py);
		ymin = Math.min(ymin, a[0]);
		ymax = Math.max(ymax, a[1]);
		if (s[1]<0) sLabel = "y="+IJ.d2s(s[0],dec)+" "+IJ.d2s(s[1],dec)+"x";
		else sLabel = "y="+IJ.d2s(s[0],dec)+" +"+IJ.d2s(s[1],dec)+"x";
		R2label = "R="+IJ.d2s(R,dec)+", "+"R2="+IJ.d2s(sqr(R),dec);
		Plot p = new Plot("Correlation plot", lab1,lab2, px, py);
		p.setLimits(xmin,xmax,ymin,ymax);

		if (marker==0) p.addPoints(a1,a2, PlotWindow.X);
		if (marker==1) p.addPoints(a1,a2, PlotWindow.CROSS);
		if (marker==2) p.addPoints(a1,a2, PlotWindow.CIRCLE);
		if (marker==3) p.addPoints(a1,a2, PlotWindow.TRIANGLE);
		if (marker==4) p.addPoints(a1,a2, PlotWindow.BOX);

		p.addLabel(0.1, 0.01, R2label);
		p.addLabel(0.7, 0.01, sLabel);
		return p.getProcessor();
	}

    public int[] s2ints(String s) {
        StringTokenizer st = new StringTokenizer(s, ", \t");
        int nInts = st.countTokens();
        int[] ints = new int[nInts];
        for(int i=0; i<nInts; i++) {
            try {ints[i] = Integer.parseInt(st.nextToken());}
            catch (NumberFormatException e) {IJ.write(""+e); return null;}
        }
        return ints;
    }

    	void writeResults(float i, float m, float s){//, int n) {

		int height=400, width=200;
		cheadings="Intervals"+"\t"+"Average"+"\t"+"Std";
		cLine =i+"\t"+IJ.d2s(m,dec)+"\t"+IJ.d2s(s,dec);
		if (ctw==null) {
			ctitle = "Increasing interval.Image correlation. Local region size = "+IJ.d2s(mSize,0)+" pixels";
			ctw = new TextWindow(ctitle,cheadings,cLine, width, height);
		}
		else
		ctw.append(cLine);

	}
    
	void writeResults(double[] l, String l1, String l2){//, int n) {

		int height=400, width=600;
		headings="Target"+"\t"+"Source"+"\t"+"R"+"\t"+"R2"+"\t"+"\t"+"Slope"+"\t"+"C"+"\t";
		aLine =l1+"\t"+l2+"\t"+IJ.d2s(l[2],dec)+"\t"+IJ.d2s(sqr(l[2]),dec)+"\t"+IJ.d2s(l[1],dec)+"\t"+IJ.d2s(l[0],dec)+"\t";
		if (tw==null) {
			title = "Image correlation. Local region size = "+IJ.d2s(mSize,0)+" pixels";
			tw = new TextWindow(title,headings,aLine, width, (height));
		}
		else
		tw.append(aLine);

	}

	public float[] getNewArray(ImageProcessor ipTemp, int mSize, int stat){
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		int nY=(int)(hh/mSize);
		int nX=(int)(ww/mSize);
		float[] array= new float[nX*nY];
		float[] localArea= new float[mSize*mSize];
		int counter=0;
 		for (int i = 0;i<nY;i++){
			int yStep=i*mSize;

			for (int ii = 0;ii<nX;ii++){
				int xStep=ii*mSize;
				float mValue = 0;
				int index =0;
				//Accumulates stats in each subimages
				for (int y = yStep; y < (yStep+mSize); y++){
					for (int x = xStep; x < (xStep+mSize); x++){
						localArea[index]=ipTemp.getPixelValue(x,y);
						index++;
					}
				}
				if (stat==4){
					FloatProcessor temp= new FloatProcessor(mSize,mSize,localArea,null);
					array[counter] = (float)getFacetDetails(temp);
				}
				else array[counter] = (float)getStatistics(localArea,mSize,stat);
				counter++;
			}
		}
		return array;
	}

	float getR(float[] d1, float[] d2){
		float t1=0, t2=0, sum=0;
		float xMean=getMean(d1);
		float yMean=getMean(d2);
		float xStd=getStd(xMean, d1);
		float yStd=getStd(yMean, d2);
		for (int i=0; i<d1.length;i++){
			t1=(d1[i]-xMean)/xStd;
			t2=(d2[i]-yMean)/yStd;
			sum=sum+(t1*t2);
		}
		float r= sum/(d1.length-1);
		return (r);
	}

	//Based on the  Facetorientation plugin by Bob Dougherty
	public float getFacetDetails(ImageProcessor ipTemp){
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		int ps=1;	//pixelSize, no calibration
		
		FloatProcessor fSlice;
		float[] gradX = (float[])getPixelValues(ipTemp);
		float[] gradY = (float[])getPixelValues(ipTemp);
		//sArea=0;
		float C=0,S=0,R=0,aT=0;
		//Put the pixel arrays into image processors for the convolutions.
		ImageProcessor ipGradX = new FloatProcessor(ww,hh);
		ipGradX.setPixels(gradX);
		ImageProcessor ipGradY = new FloatProcessor(ww,hh);
		ipGradY.setPixels(gradY);

		//Evaluate the gradients.
		ipGradX.convolve3x3(kernelX);
		ipGradY.convolve3x3(kernelY);
		float[] phi = new float[ww*hh];

		//RPD limits changed in the loops below
		for (int y = 1; y < (hh-1); y++){
			for (int x = 1; x < (ww-1); x++){
				int index = x + ww*y;
				//The derivative requires the factors 1/(2*dx) for the centered
				//finite difference and 1/4 for the coefficients of the kernel.
				double gx = gradX[index]/(8*ps);
				double gy = gradY[index]/(8*ps);
				phi[index] = (float)(Math.atan2(gy,gx)+Math.PI);
				//Calculates the mean resultant vector (R)
				//See Tovey et al., 1992
				//RPD Use double angles to avoid cancellation
				C+= Math.cos(2*phi[index]);
				S+= Math.sin(2*phi[index]);
				phi[index]=(float)(Math.toDegrees(phi[index]));if (phi[index]<0) phi[index]=360+phi[index];
			}
		}
		C=C/((hh-2)*(ww-2));S=S/((hh-2)*(ww-2));
		R=(float)(Math.sqrt(sqr2(C)+sqr2(S)));
		//RPD
		//Undo double angle									   90
		aT = (float)(Math.toDegrees(Math.atan2(S,C)/2)); 
		if (aT<0) aT=90-Math.abs(aT); 			//angles given between 0 and 180, where 180  +   0
		else aT=aT+90;
		return aT;
	}

	private float[] getPixelValues(ImageProcessor ipTemp){
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		float[] pValues = new float[ww*hh];
		for (int y = 0; y < hh; y++){
			for (int x = 0; x < ww; x++){
				int index = x + ww*y;
				pValues[index]=ipTemp.getPixelValue(x,y);
			}
		}
		return pValues;
	}
	
	public float getStatistics(float[] ipTemp, int ww, int p){
		int hh=ww;//ipTemp.getHeight();
		int N=ww*hh;
		float min = Float.MAX_VALUE;
		float max = -Float.MAX_VALUE;
		float[] rValues = new float[7];
		float ra = 0, rq = 0, sk = 0, ku = 0;
		float zMean=getMean(ipTemp);
		float zStd = getStd(zMean,ipTemp);
		int index=0;
		for (int y = 0; y < hh; y++){
			for (int x = 0; x < ww; x++){
				//float temp = ipTemp.getPixelValue(x,y);
				float temp = ipTemp[index];
				float zTemp=temp-zMean;
				sk += sqr3(zTemp/zStd);
				ku += sqr4(zTemp/zStd);
				ra += Math.abs(zTemp);
				if (temp<min) min = temp;
				if (temp>max) max = temp;
				index++;
			}
		}
		// Calculate Roughness stats (Ra, Rq, Rsk, Rku)
		rValues[0] = (float)(zMean);
		rValues[1] = (float)(zStd);
		rValues[2] = (float)(sk/N);
		rValues[3] = (float)((ku/N)-3);
		rValues[4] = min;
		rValues[5] = max;
		rValues[6] = max-Math.abs(min);
		return rValues[p];
	}

	double sqr2(double x) {return x*x;}
	double sqr3(double x) {return x*x*x;}
	double sqr4(double x) {return x*x*x*x;}
	
	double sqr(double x) {return x*x;}

	float getMean(float[] dataset){
		double mValue=0;
		for (int j=0; j<dataset.length; j++){mValue += dataset[j];}
		return (float) (mValue/dataset.length);
	}

	float getStd(float mValue,float[] dataset){
		float sValue=0;
		if (dataset.length==1) {return (float) (sValue);}
		else{
			for (int j=0; j<dataset.length; j++){sValue += sqr(mValue-dataset[j]);}
			return (float) (Math.sqrt(sValue/(dataset.length-1)));
		}
	}

	public void itemStateChanged(ItemEvent e) {
		if(e.getItemSelectable() == cbFix){
			if (e.getStateChange() == ItemEvent.SELECTED) {
				lbxFix.setVisible(true);
				tfxFix.setVisible(true);
				lbyFix.setVisible(true);
				tfyFix.setVisible(true);

			}
			else {
				lbxFix.setVisible(false);
				tfxFix.setVisible(false);
				lbyFix.setVisible(false);
				tfyFix.setVisible(false);
			}
		}
	}
}