import ij.*;
import ij.plugin.filter.*;
import ij.process.*;
import ij.gui.*;
import java.awt.*;
import java.util.*;
import ij.measure.*;
import ij.text.*;
import ij.gui.Roi.*;
import ij.process.ImageProcessor.*;
import java.awt.event.*;
import ij.io.*;
import ij.plugin.filter.RankFilters;


/*This plugin provides several routines to characterize surfaces. The plugin is based on the previous plugins for surface characterisation
Roughness_calculation, Facet_orientation, Roughness_waviness and Surface_Leveling.
The plugins can be used with any topographical image acquired with e.g.  laser profilometer, white light interferometer, AFM etc.
The input is a 32 bit image or stack in which the pixel values represent distance, z, to a surface.
The user also needs to provide the transverse spacing between the pixels, expressed in the same units as the range data.

ImageJ plugins are developed to level and assess surface roughness statistics including root mean square deviation, skewness and kurtosis.
Gaussian filtering is used to separate the waviness from the roughness and thus assess the structure at different wavelengths.
FFT bandpass filtering (Wayne Rasband and Joachim Walter) is also built in to assess the structure at different wavelengths/bands.
A routine for detecting the local facet orientation using Sobel operators is also available.

For more details see:

Chinga, G, Gregersen, Ø, Dougherty, B. "Paper surface characterisation by laser profilometry and image analysis",
J. of Microscopy and Analysis, July. 2003.
Chinga, G. and Dougherty, R.: "Quantification of surface structures", Proceedings from the ImageJ conference, Luxembourg, May 18-19 (2006).

Gary Chinga.

030708	Version 1b
031111	Version 1c
		The image labels are added to the result table.
		Updated the doFFTFiltering method.
031124	The Gaussian radius has been adjusted to half the lStruct value.
040316	Display a polar plot based on the azimuthal images (v. 1d)
050326	Roughness statistics modified, no 0-level requirement.(1e)
060120	Computes the surface area based on a suggestion given by Bob Dougherty, Optinav.com.(1f)
060204  Fixed a small bug in the calculation of the black area fraction of polar facets.(1g)
	Implemented a domain segmentation based on a gradient analysis. The domains size can be defined with
	the kernel size of a median filter.
	Implemented the calculation of the facet direction and mean resultant vector based on a paper by Curray (1956)  
	Corrected the phi value of the facets so the azimuth angle of each facet will be given according to the polar plot coordinates(v. 1kl)
060518	In Addition to the R-family, the following parameters are included, 
	FPO: Facet polar orientation relative to the plane, 
	FPOV: Standard deviation of polar orientation values
	FAD: Facet azimuthal direction. The structure orientation will be 90 degrees to the FAD values.
	MRV: Mean resultant vector. The higher the values the higher the preferred azimuthal direction. See also Curray (1956)  
	SA: Surface area in calibrated units.
	Domain segmentation in the polar (p) and azimuthal images (a). The values are given in percentage for every class
071009	SA is calculated relative to the assessed area
	Minor changes to the image label registration.
080407	Included the Rc parameter, average height of an unleveled surface. (v. 1p)	
080507	Limit the upper structure size to a 1/4 of the image size.(1q).
*/

public class SurfCharJ_1q implements PlugInFilter, Measurements,ItemListener  {

	public static final int MEDIAN=0, MEAN=1, MIN=2, MAX=3, MEANMAXMIN=4;
	ImagePlus imp;
	int i,w,h;
	double stepSize=0, pSize,lStruct,uStruct, kSize;
	double[] pLim, aLim;
	float min,max;
	float sArea,C,S,R,aT;
	boolean levelSurf=false, localSurf=false, facetAnalysis=true,medianFilter=false;
	boolean filterSurf=false, dAzimuthal=false, dPolar=false, pThresh=false, aThresh=false,entireStack=false, polarPlot= false, fPolar=false;
	boolean gFilter=false, fftFilter=false, includeInfo=true;
	boolean canceled;
	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};
	String[] pLabel = new String[15];
	float mAngle,medianAngle, sAngle;// whiteArea;
	ByteProcessor ipPlot;
	static int par=4;
	int www=400, hhh=400;
	ImageStack imsPhi, imsPolarPlot, imsVectorImage;
 	TextWindow tw;
	String title,headings,line, info;

	//Variables in input dialog
	private static String[] items = {"GF, display roughness image","GF, display waviness image", "FFT, display bandpass image"};
	protected static final int ROUGHNESS=0,WAVINESS=1,BANDPASS=2 ;
	Choice chFilters = new Choice();
	int filter=0;

	Checkbox cbLevel = new Checkbox("Level surface",levelSurf);
	Checkbox cbRoughnessAnalysis = new Checkbox("Local roughness analysis",localSurf);
	Checkbox cbFacetAnalysis = new Checkbox("Perform gradient analysis",facetAnalysis);
	Checkbox cbFilterSurf = new Checkbox("Filter surface by: ",filterSurf);
	Checkbox cbDisplayAzimuthal = new Checkbox("Display azimuthal image",dAzimuthal);
	Checkbox cbDisplayPolarPlot = new Checkbox("Display plot of azimuthal facets",polarPlot);
	Checkbox cbDisplayPolar = new Checkbox("Display polar image ",dPolar);
	Checkbox cbFilterPolar = new Checkbox("Median filter, radius (pixels)=",fPolar);

	Checkbox cbPolarThresh = new Checkbox("Threshold polar (0<t<90) ",pThresh);
	Checkbox cbAzimuthalThresh = new Checkbox("Threshold azimuthal (0<t<360) ",aThresh);

	Checkbox cbEntireStack = new Checkbox("Process entire stack ",entireStack);
	Checkbox cbIncludeInfo = new Checkbox("Include processing information ",includeInfo);

	TextField tfStep = new TextField("200",0);
	TextField tfKernel = new TextField("9",0);

	TextField tfPolarThresh = new TextField("0,15,30,45,60,75,90",0);
	TextField tfAzimuthalThresh = new TextField("0,45,90,135,180,225,270,315,360",0);

	TextField tfLowerFilter = new TextField("5",0);
	TextField tfUpperFilter = new TextField("10",0);
	TextField tfMessage = new TextField("",0);
	TextField tfMinMax = new TextField("",0);
	TextField tfMinMaxStep = new TextField("",0);
	TextField tfThreshold = new TextField("",0);
	Label lblStep = new Label("Sampling length:");
	Label lblLowerFilter = new Label("Lower structure size limit: ");
	Label lblUpperFilter = new Label("Upper structure size limit: ");
	Label lblMessage = new Label("Test");
	Label lblMinMax = new Label("Min Max");
	Label lblMinMaxStep = new Label("Min Max");
	Label lblNone = new Label("");
	Label lblNone1 = new Label("");
	Label lblNone2 = new Label("");
	Label lblNone3 = new Label("");
	Label lblNone4 = new Label("");
	Label lblNone5 = new Label("");
	Label lblNone6 = new Label("");
	Label lblNone7 = new Label("");
	Label lblNone8 = new Label("");
	Label lblNone9 = new Label("");
	Label lblNone10 = new Label("");
	Label lblPolar = new Label("");
	Label lblFPolar = new Label("");
	Label lblGradientTitle = new Label("Gradient analysis");
	Label lblGradient = new Label("");
	Label lblRosePlot = new Label("");

	public int setup(String arg, ImagePlus imp) {
		if (IJ.versionLessThan("1.37"))
			return DONE;
		this.imp = imp;
		return DOES_32;
	}

	public void run(ImageProcessor ip) {
		pLabel[0]="Rq"; pLabel[1]="Ra";pLabel[2]="Rsk";pLabel[3]="Rku";pLabel[4]="Rv";pLabel[5]="Rp";pLabel[6]="Rt";pLabel[7]="Rc";
		pLabel[8]="FPO"; pLabel[9]="MFOV";pLabel[10]="FAD"; pLabel[11]="MRV";pLabel[12]="SA"; 
		String units = imp.getCalibration().getUnits();
		float temp = 0;
		int parameters=8;
		int bPar=8;
		ImageStack stack = imp.getStack();
		int w = ip.getWidth();
		int h = ip.getHeight();
		int c;

		if (h<w) c=h; else c=w;		//Calculates the smallest dimension
		int nSlices = imp.getStackSize();
		int nMeasurements=0;

		Calibration cal = imp.getCalibration();
		pSize = cal.pixelWidth;

		getDetails(c,pSize, nSlices);
		if (canceled) return;
		//Sums up the number of variables necessary in the result array.
		if (facetAnalysis){
			parameters += 5;		
			if (pThresh) parameters +=pLim.length-1;
			if (aThresh) parameters +=aLim.length-1;
		}
		
		if (entireStack) nMeasurements=nSlices;
		else nMeasurements=1;
		String[] imageLabel = new String[nMeasurements];

		int ms= (int)(stepSize/pSize);
		lStruct=lStruct/pSize;
   		uStruct=uStruct/pSize;
   		registerInfo(units);
		//array for results
		float[][] roughnessValues= new float[nSlices][parameters];
		float[] tRValues=new float[bPar];
		float[] tFractions=new float[pLim.length-1];
		float[] aFractions=new float[aLim.length-1];

		ImageStack imsTheta = new ImageStack(w,h);
		ImageStack imsPhi = new ImageStack(w,h);
		ImageStack imsPolarPlot = new ImageStack(www,hhh);
		ImageStack imsThreshold = new ImageStack(w,h);
		ImageStack imsAziThreshold = new ImageStack(w,h);

		FloatProcessor ipThresh=new FloatProcessor(w,h);	//contains the polar thresholded image
		FloatProcessor ipAThresh=new FloatProcessor(w,h);	//contains the azimuthal thresholded image
		ImageProcessor ipSlice;//, ipFiltSlice;

		for (int n = 0; n < nMeasurements; n++){
			int basicPar=bPar;
			if (facetAnalysis) basicPar+=5;
			IJ.showProgress((double)(n+1)/nSlices);
			IJ.showStatus("a: "+n+"/"+nSlices);
			//imageLabel[n]= stack.getSliceLabel(n+1);
			imageLabel[n]= stack.getShortSliceLabel(n+1);
			if (entireStack){
				imp.setSlice(n+1);
				ipSlice = stack.getProcessor(n+1);
			}
			else ipSlice = stack.getProcessor(imp.getCurrentSlice());
			ipSlice.convertToFloat();

			if (levelSurf) ipSlice = levelSurface(ipSlice);

			if (filterSurf){
					if (filter==ROUGHNESS)
						ipSlice=doGaussianFiltering(ipSlice, (lStruct/2),true);		//Gaussian radius is equal to half the lStruct size.
					else if (filter==WAVINESS)
						ipSlice=doGaussianFiltering(ipSlice, (lStruct/2),false);
					else if (filter==BANDPASS)
						ipSlice=doFFTFiltering(ipSlice,lStruct,uStruct);
			}

			if (!localSurf) tRValues = getRoughnessValues(ipSlice);
			else tRValues = getSubRoughnessValues(ipSlice,ms);
			for (int i=0;i<bPar;i++)
				roughnessValues[n][i]=tRValues[i];		

			if (facetAnalysis){
				FloatProcessor ipPolar=getFacetDetails(ipSlice,pSize,true);
				FloatProcessor ipAzimuthal = getFacetDetails(ipSlice,pSize,false);
				if(pThresh) {
					ipThresh=thresholdPolar(ipPolar,pLim,medianFilter,kSize);
					tFractions=getAreaFractions(ipThresh,pLim,kSize);
					for (int i=0;i<tFractions.length;i++)
						roughnessValues[n][basicPar+i]=tFractions[i];
					basicPar +=pLim.length-1;
				}
				if(aThresh) {
					ipAThresh=thresholdPolar(ipAzimuthal,aLim,medianFilter,kSize);
					aFractions=getAreaFractions(ipAThresh,aLim,kSize);
					for (int i=0;i<aFractions.length;i++)
						roughnessValues[n][basicPar+i]=aFractions[i];
				}
				roughnessValues[n][8]= calculateMean(ipPolar);
				roughnessValues[n][9]= calculateStd(ipPolar,roughnessValues[n][7]);
				roughnessValues[n][10]= aT;
				roughnessValues[n][11]= R;
				roughnessValues[n][12]=sArea; //surface area
				if (dPolar||pThresh){
					ipPolar.setMinAndMax(0,0);
					imsTheta.addSlice("Polar  image (Theta)",ipPolar);
				}
				if (pThresh) imsThreshold.addSlice("",ipThresh);

				if (dAzimuthal){
					ipAzimuthal.setMinAndMax(0,0);
					imsPhi.addSlice("Azimuthal  image of "+imageLabel[n],ipAzimuthal);
				}
				if (aThresh) imsAziThreshold.addSlice("",ipAThresh);
				if (polarPlot){
					float[] phi = (float[])ipAzimuthal.getPixelsCopy();
					int[] aAngles=registerPolarAngles(phi,w,h);		//register the frequencies of angles 0-360
					int[][] xyCoords=calCulatePolarPlot(aAngles, www, hhh);			//returns the x,y coordinates in a polar plot
					ipPlot=makePolarImage(xyCoords,www,hhh);
					imsPolarPlot.addSlice(imageLabel[n],ipPlot);
				}

			}
			ipSlice.setMinAndMax(0,0);
			imp.setProcessor(imp.getShortTitle(), ipSlice);
		}
		if(dAzimuthal) createImagePlus(imsPhi, "Azimuthal images", cal, kSize, w,h);
		if(dPolar) createImagePlus(imsTheta,"Polar images", cal,kSize, w,h);
		if(polarPlot) createImagePlus(imsPolarPlot, "Polar plot images", cal, kSize, w,h);
		if(pThresh) createImagePlus(imsThreshold,"Polar coded images", cal, kSize, w,h);
		if(aThresh) createImagePlus(imsAziThreshold,"Azimuthal coded images", cal, kSize, w,h);

		imp.updateAndDraw();
		if (filter==ROUGHNESS||filter==WAVINESS) IJ.run("Fire");
		if (includeInfo) showInfo(info);
		writeResults(roughnessValues, nMeasurements, pLabel,parameters,title, imageLabel,pLim,aLim);
		IJ.showProgress(1.0);
	}

	void getDetails(int mMask, double ps, int numSlices) throws NumberFormatException{
		GenericDialog gd = new GenericDialog("Surface characterisation 1kl...");
		String units = imp.getCalibration().getUnits();
		gd.setLayout(new GridLayout(16,1));
		gd.add(cbLevel);gd.add(lblNone);
		gd.add(cbRoughnessAnalysis);gd.add(lblNone1);
		gd.add(lblStep);gd.add(tfStep);
		gd.add(lblMinMaxStep);gd.add(tfMinMaxStep);
		gd.add(cbFacetAnalysis);
		gd.add(lblGradient);
		gd.add(cbDisplayPolar);gd.add(cbDisplayAzimuthal);
		gd.add(cbFilterPolar);gd.add(tfKernel);

		gd.add(cbPolarThresh);gd.add(tfPolarThresh);
		gd.add(cbAzimuthalThresh);gd.add(tfAzimuthalThresh);
		gd.add(cbDisplayPolarPlot);gd.add(lblRosePlot);
		gd.add(cbFilterSurf);
		chFilters.add(items[0]);
		chFilters.add(items[1]);
		chFilters.add(items[2]);
		gd.add(chFilters);
		gd.add(lblLowerFilter);gd.add(tfLowerFilter);
		gd.add(lblUpperFilter);gd.add(tfUpperFilter);
		gd.add(lblMinMax);gd.add(tfMinMax);
		gd.add(cbIncludeInfo);//gd.add(lblNone8);
		gd.add(cbEntireStack);//gd.add(lblNone7);
		//gd.add(lblMessage);gd.add(tfMessage);

		lblMinMaxStep.setVisible(false);tfMinMaxStep.setVisible(false);
		lblMinMax.setVisible(false);tfMinMax.setVisible(false);
		lblStep.setVisible(false);tfStep.setVisible(false);
		//lblMessage.setVisible(false);tfMessage.setVisible(false);
		lblLowerFilter.setVisible(false);tfLowerFilter.setVisible(false);
		lblUpperFilter.setVisible(false);tfUpperFilter.setVisible(false);
		chFilters.setVisible(false);
		cbDisplayAzimuthal.setVisible(true);cbDisplayPolarPlot.setVisible(true);
		cbDisplayPolar.setVisible(true);
		cbPolarThresh.setVisible(true);tfPolarThresh.setVisible(true);
		cbFilterPolar.setVisible(true);tfKernel.setVisible(true);
		cbAzimuthalThresh.setVisible(true);tfAzimuthalThresh.setVisible(true);
		
		tfThreshold.setVisible(false);
		if (numSlices>1) cbEntireStack.setVisible(true);
		else cbEntireStack.setVisible(false);
		lblStep.setText("Sampling length ("+units+"): ");
		lblMinMaxStep.setText("Min="+(int)(4*ps)+", Max="+(int)(mMask*ps/2)+units);
		lblLowerFilter.setText("Lower structure size limit ("+units+") : ");
		lblUpperFilter.setText("Upper structure size limit ("+units+") : ");
		lblMinMax.setText("Max="+(int)(mMask*ps/4)+units);
		tfMinMax.setVisible(false);

		cbRoughnessAnalysis.addItemListener(this);
		cbPolarThresh.addItemListener(this);
		cbFilterSurf.addItemListener(this);
		gd.add(lblNone3);
		gd.showDialog();

		if (gd.wasCanceled()) {
			canceled = true;
			return;
		}

		levelSurf = cbLevel.getState();
		localSurf = cbRoughnessAnalysis.getState();
		facetAnalysis = cbFacetAnalysis.getState();
		if (facetAnalysis){
			dAzimuthal = cbDisplayAzimuthal.getState();
			polarPlot= cbDisplayPolarPlot.getState();
			dPolar = cbDisplayPolar.getState();
			pThresh = cbPolarThresh.getState();
			aThresh = cbAzimuthalThresh.getState();
			medianFilter = cbFilterPolar.getState();
		}

		filterSurf = cbFilterSurf.getState();
		includeInfo=cbIncludeInfo.getState();
		entireStack = cbEntireStack.getState();
		filter = chFilters.getSelectedIndex();

		if(localSurf){
			try{
				stepSize = (double)(Double.parseDouble(tfStep.getText()));

			} catch(NumberFormatException e){
				IJ.showMessage("Invalid input number, R-values will be based on whole image");
				localSurf = false;
			}
		}

		if(pThresh||aThresh){
			try{
				kSize = (double)(Double.parseDouble(tfKernel.getText()));
				if (!medianFilter) kSize=0;		//0 for not affecting the image dimensions
				if (pThresh) pLim = s2doubles(tfPolarThresh.getText()); else pLim = s2doubles("0");
				if (aThresh) aLim = s2doubles(tfAzimuthalThresh.getText()); else aLim = s2doubles("0");

			} catch(NumberFormatException e){
				IJ.showMessage("Invalid input number, thresholding will not be performed");
				pThresh = false;
				aThresh = false;
			}
		}
		else {
			pLim = s2doubles("0");
			aLim = s2doubles("0");
		}

		if(filterSurf){
			try{
				lStruct = (double)(Double.parseDouble(tfLowerFilter.getText()));
				uStruct = (double)(Double.parseDouble(tfUpperFilter.getText()));
				if ((uStruct<lStruct)&&(filter==BANDPASS)){ 								//verify if the input values are valid
					uStruct=lStruct*2;
					IJ.showMessage("Upper size limit is set to "+uStruct+" when FFT filtering");
				}
				if (uStruct>(int)(mMask*ps)/4){
					IJ.showMessage("Invalid input number, Filtering will not be performed");
					filterSurf = false;
				}
				if (filter==ROUGHNESS||filter==WAVINESS){
					if (lStruct>(int)(mMask*ps)/4){
						IJ.showMessage("Invalid input number, Filtering will not be performed");
						filterSurf = false;
					}
				}

			} catch(NumberFormatException e){
				IJ.showMessage("Invalid input number, Filtering will not be performed");
				filterSurf = false;
			}
		}

		if (localSurf && (stepSize<(int)(4*ps)) || (stepSize>(mMask*ps)/2)) {
			canceled = true;
			IJ.showMessage("Min and Max recommended mask size are " +(int)(4*ps)+" and "+(int)(mMask*ps/2)+" respectively");
			return;
		}

	}

	public void itemStateChanged(ItemEvent e) {
		if(e.getItemSelectable() == cbRoughnessAnalysis){
			if (e.getStateChange() == ItemEvent.SELECTED) {
				lblStep.setVisible(true);
				tfStep.setVisible(true);
				lblMinMaxStep.setVisible(true);
				tfMinMaxStep.setVisible(false);
			}
			else {
				lblStep.setVisible(false);
				tfStep.setVisible(false);
				lblMinMaxStep.setVisible(false);
				tfMinMaxStep.setVisible(false);
			}
		}

		else{
		
					if(e.getItemSelectable() == cbFilterSurf){
						if (e.getStateChange() == ItemEvent.SELECTED) {
							lblMinMax.setVisible(true);
							lblLowerFilter.setVisible(true);
							tfLowerFilter.setVisible(true);
							lblUpperFilter.setVisible(true);
							tfUpperFilter.setVisible(true);
							chFilters.setVisible(true);
						}
						else {
							lblMinMax.setVisible(false);
							lblLowerFilter.setVisible(false);
							tfLowerFilter.setVisible(false);
							lblUpperFilter.setVisible(false);
							tfUpperFilter.setVisible(false);
							chFilters.setVisible(false);
						}
					}
		}
	}

	void registerInfo(String u){
		info="Processing information: ";
		if (localSurf)info+="Sampling length="+stepSize+u;
		else info+="No sampling";
		if (levelSurf) info+=", Surface leveling";
		if (filterSurf){
				if (filter==ROUGHNESS) info+=", Gaussian roughness, r="+(lStruct/2*pSize)+u;
				if (filter==WAVINESS) info+=", Gassian waviness, r="+(lStruct/2*pSize)+u;
				if (filter==BANDPASS) info+=", FFT bandpass ("+(lStruct*pSize)+" and "+(uStruct*pSize)+u+")";
		}
		if (medianFilter) info+=", Median filter, radius= "+kSize;
	}

	FloatProcessor doFFTFiltering(ImageProcessor ipSlice, double lS, double uS){
		ImagePlus tempImage= new ImagePlus("Temp image", ipSlice);
		WindowManager.setTempCurrentImage(tempImage);
		IJ.run("Bandpass Filter...", "filter_large="+uS+" filter_small="+lS+" suppress=None tolerance=2");
		FloatProcessor ipTemp=(FloatProcessor)tempImage.getProcessor();
		return ipTemp;
	}

	FloatProcessor doGaussianFiltering(ImageProcessor ipTemp, double r, boolean roughness){
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		float[] wArray = (float[])getPixelValues(ipTemp);
		FloatProcessor ipWaviness= new FloatProcessor(ww,hh);
		ipWaviness.setPixels(wArray);

		float[] rArray = (float[])getPixelValues(ipTemp);
		FloatProcessor ipRoughness= new FloatProcessor(ww,hh);
		ipRoughness.setPixels(rArray);
		GaussianBlur gb = new GaussianBlur();
		gb.blur(ipWaviness,r);
		if (roughness){
			for (int y = 0; y < hh; y++){
				for (int x = 0; x < ww; x++){
					int index = x + ww*y;
						rArray[index] = rArray[index]-wArray[index];
				}
			}
		}
		if (roughness) return ipRoughness;
		else return ipWaviness;
	}


	float[] getSubImage(float[] tPhi, int yIndex, int xIndex,int subMask,int imageWidth){
		float[] tempPhi = new float[subMask*subMask];
		for (int yy=yIndex; yy<(subMask+yIndex); yy++){
			for (int xx=xIndex; xx<(subMask+xIndex);xx++){
				int pIndex = xx + imageWidth*yy;
				int tIndex = ((xx-xIndex) + subMask*(yy-yIndex));
				tempPhi[tIndex]= tPhi[pIndex];
			}
		}
		return tempPhi;
	}

	public float[] getRoughnessValues(ImageProcessor ipTemp){
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		int N = ww*hh;

		float min = Float.MAX_VALUE;
		float max = -Float.MAX_VALUE;
		float[] rValues = new float[8];
		float ra = 0, rq = 0, sk = 0, ku = 0, rc=0;
		float zMean=calculateMean(ipTemp);
		float zStd = calculateStd(ipTemp,zMean);

		for (int y = 0; y < hh; y++){
			for (int x = 0; x < ww; x++){
				float temp = ipTemp.getPixelValue(x,y);
				float zTemp=temp-zMean;
				sk += sqr3(zTemp/zStd);
				ku += sqr4(zTemp/zStd);
				ra += Math.abs(zTemp);
				rc += temp;
				if (temp<min) min = temp;
				if (temp>max) max = temp;
			}
		}
		// Calculate Roughness stats (Ra, Rq, Rsk, Rku)
		rValues[1] = (float)(ra/(N));
		rValues[0] = (float)(zStd);
		//float Pq= rValues[0];

		rValues[2] = (float)(sk/N);
		rValues[3] = (float)((ku/N)-3);
		rValues[4] = min;
		rValues[5] = max;
		rValues[6] = max-min;
		rValues[7] = (float)(rc/(N));
		
		return rValues;
	}

	public float[] getSubRoughnessValues(ImageProcessor ipTemp, int mSize){
		//counter=0;
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		int N=(int)sqr2(mSize);
		int nY=(int)(hh/mSize);
		int nX=(int)(ww/mSize);
		float[] rValues= new float[8];
		float[] tempValues = new float[mSize*mSize];

		for (int i = 0;i<nY;i++){
			int yStep=i*mSize;
			for (int ii = 0;ii<nX;ii++){
				int xStep=ii*mSize;
				float min = Float.MAX_VALUE;
				float max = -Float.MAX_VALUE;
				float ra = 0, rq = 0, sk = 0, ku = 0, rc=0;
				int counter=0;
				for (int y = yStep; y < (yStep+mSize); y++){
					for (int x = xStep; x < (xStep+mSize); x++){
						tempValues[counter]=ipTemp.getPixelValue(x,y);
						counter++;
					}
				}
				FloatProcessor floatTemp = new FloatProcessor(mSize,mSize,tempValues,null);
				float zMean=calculateMean(floatTemp);
				float zStd = calculateStd(floatTemp,zMean);
				//Accumulates stats in each subimages
				for (int y = yStep; y < (yStep+mSize); y++){
					for (int x = xStep; x < (xStep+mSize); x++){
						float temp = ipTemp.getPixelValue(x,y);
						float zTemp=temp-zMean;
						sk += sqr3(zTemp/zStd);
						ku += sqr4(zTemp/zStd);
						ra += Math.abs(zTemp);
						rc += temp;
						
						if (temp<min) min = temp;
						if (temp>max) max = temp;
					}
				}
				// Calculate Roughness stats of the subimage(Ra, Rq, Rsk, Rku)
				float Pa = (float)(ra/N);
				float Pq = zStd;
				float Psk = (float)(sk /N);
				float Pku = (float)((ku / N)-3);
				float Pc = (float) (rc/N);

				//sums the statistics of each subimage
				rValues[0] += Pq;
				rValues[1] += Pa;
				rValues[2] += Psk;
				rValues[3] += Pku;
				rValues[4] += min;
				rValues[5] += max;
				rValues[7] += Pc;
				
			}
		}

		//Calculates the average of all the subimages for each image
		rValues[0] = rValues[0]/(nX*nY);
		rValues[1] = rValues[1]/(nX*nY);
		rValues[2]=rValues[2]/(nX*nY);
		rValues[3] = rValues[3]/(nX*nY);
		//Assign the values of min and max to the variables Rp and Rv
		rValues[4] = rValues[4]/(nX*nY);
		rValues[5] = rValues[5]/(nX*nY);
		rValues[6] = rValues[6] + (float)(rValues[5] - rValues[4]);
		rValues[7] = rValues[7]/(nX*nY);

		return rValues;
	}

	public FloatProcessor levelSurface(ImageProcessor ipTemp){
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		int n;
		float Ex1, Ex1x1, Ex2x2, Ex1x2, Ey, Ex2;
		float x1Mean, x2Mean, yMean, Eyy, Ex1y, Ex2y;
		float Sx1y, Sx2y, Sx1x2, Sx1x1, Sx2x2, b1, b2, b0;

 		float[] tempArray = (float[])getPixelValues(ipTemp);
		//Runs the same routines 4 times
		for (int it=0; it<4;it++){
			n=0;
			Ex1=0; Ex1x1=0; Ex2x2=0; Ex1x2=0; Ey=0; Ex2=0;
			x1Mean = 0; x2Mean=0; yMean = 0; Eyy = 0; Ex1y = 0; Ex2y=0;
			Sx1y = 0; Sx2y=0; Sx1x2 = 0; Sx1x1=0; Sx2x2=0; b1 = 0; b2 = 0; b0 = 0;

			//Go through the pixel values in each slice (image)
			//x1 og x2 are the pixel coordinates, y are the topography values
			for (int x2= 0; x2 < hh; x2++){
				for (int x1 = 0; x1 < ww; x1++){
					n++;
					int index = x1 + ww*x2;
					//int index2 = x2 + hh*x1;
					float y = tempArray[index];

					Ex1 = Ex1 + x1;			//Sums x, y...
					Ex2 = Ex2 + x2;
					Ey = Ey + y;
					Ex1x1 = Ex1x1 + x1*x1;	//Sums of squares and
					Ex2x2 = Ex2x2 + x2*x2;
					Ex1x2 = Ex1x2 + x1*x2;	//cross products of variables
					Ex1y = Ex1y + x1*y;
					Ex2y = Ex2y + x2*y;
				}
			}

			yMean = Ey / n;
			x1Mean = Ex1 / n;
			x2Mean = Ex2 / n;
			Sx1x2 = Ex1x2 - (n*(x1Mean)*(x2Mean));
			Sx1x1 = Ex1x1 - (n*(x1Mean)*(x1Mean));
			Sx2x2 = Ex2x2 - (n*(x2Mean)*(x2Mean));
			Sx1y = Ex1y - (n*(x1Mean)*(yMean));
			Sx2y = Ex2y - (n*(x2Mean)*(yMean));

			//calculate the values of the least square estimates
			b1 = (((Sx2x2*Sx1y)+(Sx1x2*Sx2y))/((Sx1x1*Sx2x2)+(Sx1x2*Sx1x2)));
			b2 = ((Sx2y-(b1*Sx1x2))/(Sx2x2));
			b0 = yMean - b1*x1Mean - b2*x2Mean;

			for (int x2 = 0; x2 < hh; x2++){
				for (int x1 = 0; x1 < ww; x1++){
					int index = x1 + ww*x2;
					float y = tempArray[index];
					//regression model applied to each pixel in the image
					tempArray[index] = (float) (y-(b0+b1*x1 + b2*x2));
				}
			}
		}
		FloatProcessor fSlice = new FloatProcessor(ww, hh, tempArray,null);
		return fSlice;
	}

	//Based on the  Facetorientation plugin by Bob Dougherty
	public FloatProcessor getFacetDetails(ImageProcessor ipTemp, double ps, boolean polarSlice){
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		FloatProcessor fSlice;
		float[] gradX = (float[])getPixelValues(ipTemp);
		float[] gradY = (float[])getPixelValues(ipTemp);
		sArea=0;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[] theta = new float[ww*hh];
		float[] phi = new float[ww*hh];
		float[] threshold = 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);
				theta[index] = (float)(Math.atan(Math.sqrt(gx*gx + gy*gy)));
				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]);
				//Computes the surface element by da = sqrt(1 + (tan(theta))^2)
				//Suggestion given by Bob Dougherty, Optinav.com
				sArea += (float)(sqr2(ps)*Math.sqrt(1 + sqr2(Math.tan(theta[index]))));
				phi[index]=(float)(Math.toDegrees(phi[index]));if (phi[index]<0) phi[index]=360+phi[index];
				theta[index]=(float)(Math.toDegrees(theta[index]));
			}
		}
		sArea=(float)(sArea/((ww-2)*(hh-2)*sqr2(ps)));
		C=C/((hh-2)*(ww-2));S=S/((hh-2)*(ww-2));
		R=(float)(Math.sqrt(sqr2(C)+sqr2(S)));
		//RPD
		//Undo double angle
		aT = (float)(Math.toDegrees(Math.atan2(S,C)/2));
		if (aT<0) aT=180+aT; //orientation between 0 and 360.
		if (polarSlice)  fSlice = new FloatProcessor(ww, hh, theta,null);
		else fSlice = new FloatProcessor(ww, hh, phi,null);
		return fSlice;
	}

	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;
	}

/**************************** FACET ORIENTATION - POLAR PLOT **************************/
	ByteProcessor makePolarImage(int[][] xy, int ww, int hh){
		ByteProcessor ipP = new ByteProcessor(ww,hh);
		//ipP.setColor(Color.white);
		ipP.setColor(new Color(255,255,255));ipP.fill();
		ipP.setColor(new Color(0,0,0));
		ipP.setLineWidth(1);
		ipP.setFont(new Font("SansSerif",Font.PLAIN,12));
		ipP.setAntialiasedText(false);
		ipP.moveTo(xy[0][0], xy[1][0]);
		for (int ii=0; ii<360; ii++)
			//ipP.putPixelValue(xy[0][ii], xy[1][ii], 255);
			ipP.lineTo(xy[0][ii], xy[1][ii]);

		ipP.lineTo(xy[0][0], xy[1][0]);
		//ipP.rotate(-90);						//rotattes the image since the facets are pointing 90 degrees t
											//to the structure orientation
		//Write coordinates and cross in the middle.
		ipP.moveTo(ww/2-5, ww/2);
		ipP.lineTo(ww/2+5, ww/2);
		ipP.moveTo(ww/2, ww/2-5);
		ipP.lineTo(ww/2, ww/2+5);
		ipP.drawString("0", ww-14, ww/2);
		ipP.drawString("90", ww/2-6, 13);
		ipP.drawString("180", 2, ww/2);
		ipP.drawString("270", ww/2, ww);
		ipP.setLineWidth(1);
		ipP.drawRect(ww-30, ww/2-14, 30, 14);
		ipP.drawRect(ww/2-15, 0, 40, 14);
		ipP.drawRect(0, ww/2-14, 30, 14);
		ipP.drawRect(ww/2-10, ww-14, 40, 14);
		//ipP.invert();

		return ipP;
	}

	int[] registerPolarAngles(float[] azi, int iw, int ih){
		int[] aziAngle= new int[361];
		for (int y = 1; y < ih-2; y++){
			for (int x = 1; x < iw-2; x++){
				int index = x + iw*y;
				if (azi[index]<0) azi[index]= 360+azi[index];

				//IJ.write(IJ.d2s(azi[index],2));
				aziAngle[(int)azi[index]]++;
			}
		}
		//Filters the values in order to remove deviating peaks at 0,45, 90, 135 etc...
		float[] tempAngles = new float[3];
		FloatProcessor tAngle;
		for (int ii=0;ii<358;ii++){
			for (int iii=0; iii<3; iii++)
				tempAngles[iii]=aziAngle[ii+iii];
			tAngle = new FloatProcessor(3, 1, tempAngles,null);
			aziAngle[ii]= (int)findMedian(tAngle);
		}
		tempAngles[0]=aziAngle[358]; tempAngles[1]=aziAngle[359]; tempAngles[2]=aziAngle[0];
		tAngle = new FloatProcessor(3, 1, tempAngles,null);
		aziAngle[358]=(int)findMedian(tAngle);
		tempAngles[0]=aziAngle[359]; tempAngles[1]=aziAngle[0]; tempAngles[2]=aziAngle[1];
		tAngle = new FloatProcessor(3, 1, tempAngles,null);
		aziAngle[359]=(int)findMedian(tAngle);

		return aziAngle;
	}

	int[][] calCulatePolarPlot(int[] aAngles, int ww, int hh){
		int [][] Coords = new int[2][360];
		//int [] yCoord = new int[360];
		for (int ii=0; ii<360; ii++){
			double angle = ii*Math.PI/180; //convert the angles to radians

			Coords[0][ii]= -(int)(aAngles[ii]*Math.cos(angle)); //Inverse the x-axis
			Coords[1][ii]= (int)(aAngles[ii]*Math.sin(angle)); //Inverse the x-axis

		}

		//move the plot to positive values and
		int yMin = 999999,xMin = 999999,yMax = -999999,xMax = -999999;
		for (int ii=0; ii<360;ii++){
			if (Coords[0][ii]<xMin) xMin = Coords[0][ii];
			else if (Coords[0][ii]>xMax) xMax = Coords[0][ii];
			if (Coords[1][ii]<yMin) yMin = Coords[1][ii];
			else if (Coords[1][ii]>yMax) yMax = Coords[1][ii];
		}

		//rescale to fit the 450x450 image
		double scaleFactor=1;
		double maxDist;

		if ((xMax-xMin)>(yMax-yMin)) maxDist = (xMax-xMin);
		else maxDist = (yMax-yMin);

		scaleFactor= maxDist/(ww-100);

		int xP = (int)((ww-((xMax-xMin)/scaleFactor))/2); 	//Finds the starting point
		int yP = (int)((ww-((yMax-yMin)/scaleFactor))/2);	//The plot is placed on the center of the image.

		//moves the plot to the center of the image.
		for (int ii=0; ii<360;ii++){
			Coords[0][ii] = (int)(xP+(Coords[0][ii]+ Math.abs(xMin))/scaleFactor);
			Coords[1][ii] = (int)(yP+(Coords[1][ii]+ Math.abs(yMin))/scaleFactor);
		}
		return Coords;
	}

	int[][] makePlotLabel(int[] aAngles, int ww, int hh){
		int [][] Coords = new int[2][360];
		//int [] yCoord = new int[360];
		for (int ii=0; ii<360; ii++){
			double angle = ii*Math.PI/180; //convert the angles to radians

			Coords[0][ii]= -(int)(aAngles[ii]*Math.cos(angle)); //Inverse the x-axis
			Coords[1][ii]= (int)(aAngles[ii]*Math.sin(angle)); //Inverse the x-axis

		}

		//move the plot to positive values and
		int yMin = 999999,xMin = 999999,yMax = -999999,xMax = -999999;
		for (int ii=0; ii<360;ii++){
			if (Coords[0][ii]<xMin) xMin = Coords[0][ii];
			else if (Coords[0][ii]>xMax) xMax = Coords[0][ii];
			if (Coords[1][ii]<yMin) yMin = Coords[1][ii];
			else if (Coords[1][ii]>yMax) yMax = Coords[1][ii];
		}

		//rescale to fit the 450x450 image
		double scaleFactor=1;
		double maxDist;

		if ((xMax-xMin)>(yMax-yMin)) maxDist = (xMax-xMin);
		else maxDist = (yMax-yMin);

		scaleFactor= maxDist/(ww-100);

		int xP = (int)((ww-((xMax-xMin)/scaleFactor))/2); 	//Finds the starting point
		int yP = (int)((ww-((yMax-yMin)/scaleFactor))/2);	//The plot is placed on the center of the image.

		//moves the plot to the center of the image.
		for (int ii=0; ii<360;ii++){
			Coords[0][ii] = (int)(xP+(Coords[0][ii]+ Math.abs(xMin))/scaleFactor);
			Coords[1][ii] = (int)(yP+(Coords[1][ii]+ Math.abs(yMin))/scaleFactor);
		}
		return Coords;
	}


	/**************************************************************************************************************/

	void createImagePlus(ImageStack imsTemp, String txt, Calibration c, double kSize, int ww, int hh){
		// Create new images using the new stacks.
		ImagePlus impTemp = new ImagePlus(txt,imsTemp);
		impTemp.setCalibration(c);
		impTemp.setStack(null,imsTemp);
		impTemp.show();
		WindowManager.setTempCurrentImage(impTemp);IJ.run("Fire");
		IJ.makeRectangle((int)kSize, (int)kSize, (int)(ww-2*kSize), (int)(hh-2*kSize));
	}

	//Performs domain segmentation
 	FloatProcessor thresholdPolar(ImageProcessor ipTemp,double[] t, boolean mFilter,double kSize){
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		double zero = t[0];
		RankFilters rf = new RankFilters();
		float[] tArray=new float[ww*hh];
		FloatProcessor fSlice = new FloatProcessor(ww,hh, tArray,null);

		fSlice.setPixels(ipTemp.getPixelsCopy());
		rf.rank(fSlice,kSize, MEDIAN);

		float[] floatArray=(float[])getPixelValues(fSlice);
		for (int i=0; i<t.length-1;i++){
			for (int y = 0; y < hh; y++){
				for (int x = 0; x < ww; x++){
					int index = x + ww*y;
					if (floatArray[index]==zero) 
						tArray[index]=(float)(t[1]);	
					else if (floatArray[index]>t[i] && floatArray[index]<=t[i+1]) 
						tArray[index]=(float)(t[i+1]);
				}
			}
		}
		fSlice.setPixels(tArray);
		return fSlice;
	}

	//Quantify the area fraction of each domain	
	float[] getAreaFractions(FloatProcessor ipTemp, double[] pLim,double kSize){
		float[] tFractions = new float[pLim.length-1];
		int ww=(int)(ipTemp.getWidth()-2*kSize);
		int hh=(int)(ipTemp.getHeight()-2*kSize);
		double zero = pLim[0];
		
		// Dont take into account the edges of the images due to erroneous gradients			
		for (int y = (int)kSize+1; y < (int)(hh+kSize-1); y++){		//gets area fraction within the ROI (kSize,kSize,ww-2kSize,hh-2kSize)
			for (int x = (int)kSize+1; x < (int)(ww+kSize-1); x++){
				if (ipTemp.getPixelValue(x,y)==zero)
					tFractions[0]++;
				else{
					for (int ii=0;ii<pLim.length-1;ii++)
						if (ipTemp.getPixelValue(x,y)>pLim[ii] && ipTemp.getPixelValue(x,y)<=pLim[ii+1]) tFractions[ii]++; //&&ipTemp.getPixelValue(x,y)<(pLim[ii]+1)) 				
				}
			}
		}
		for (int ii=0;ii<tFractions.length;ii++){
			tFractions[ii]=tFractions[ii]/((ww-2)*(hh-2))*100; //substract 2 due to the erroneous edge-gradients
		}
		return tFractions;
	}

	float calculateMean(ImageProcessor ipTemp){
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		float[] tempSlice = (float[])getPixelValues(ipTemp);
		double mValue=0;
		int counter=(ww*hh);
		for (int j=0; j<counter; j++){mValue += tempSlice[j];}
		return (float)(mValue/counter);
	}

	float calculateStd(ImageProcessor ipTemp,float mValue){
		int ww=ipTemp.getWidth();
		int hh=ipTemp.getHeight();
		float[] tempSlice = (float[])getPixelValues(ipTemp);
		int counter=ww*hh;
		float sValue=0;
		if (counter==1) {return (float) (sValue);}
		else{
			for (int j=0; j<counter; j++){sValue += sqr2(mValue-tempSlice[j]);}
			return (float) (Math.sqrt(sValue/(counter-1)));
		}
	}

	public double[] s2doubles(String s) {
		StringTokenizer st = new StringTokenizer(s, ", \t");
		int nDoubles = st.countTokens();
		double[] doubles = new double[nDoubles];
		for(int i=0; i<nDoubles; i++) {
			try {doubles[i] = Double.parseDouble(st.nextToken());}
			catch (NumberFormatException e) {IJ.write(""+e); return null;}
        	}
		Arrays.sort(doubles);
		return doubles;
	}



	float findMedian(ImageProcessor ipTemp){//float mValue){
		float[] input = (float[])getPixelValues(ipTemp);
		float[] input2 = new float[input.length];
		int lengthArray = input.length;
		int mid;
		float med;
		for(int k=0;k<input.length;k++){
			input2[k] = input[k];
		}
		Arrays.sort(input2);
		if((lengthArray%2)!=0){
			mid=lengthArray/2;
			med = input2[mid];
		}
		else{
			mid = (int)((lengthArray/2)+0.5);
			med = (input2[mid-1] + input2[mid])/2;
		}
		return med;
	}

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

	void writeResults(float[][] res,int nS, String [] parLabels, int par, String t, String[] iLabel,double[] pLim, double[] aLim){//, int n) {
		int height=400, width=800;
		headings="Slice"+"\t";
		//IJ.showMessage(IJ.d2s(par)+"  "+IJ.d2s(pLim.length-1)+"  "+IJ.d2s(aLim.length-1));
		for (int x=0;x<(par-(pLim.length-1)-(aLim.length-1));x++)
			headings += parLabels[x]+"\t";

		//Enters the polar fraction areas.
		//if (pThresh)
		for (int xx=0;xx<pLim.length-1;xx++)
			headings += IJ.d2s(pLim[xx],0)+":"+IJ.d2s(pLim[xx+1],0)+"p\t";
		//if (aThresh)
		for (int xx=0;xx<aLim.length-1;xx++)
			headings += IJ.d2s(aLim[xx],0)+":"+IJ.d2s(aLim[xx+1],0)+"a\t";

		for (int y=0;y<nS;y++){  //
			if (iLabel[y]==null) line =IJ.d2s((y+1),0)+"\t";
			else line = iLabel[y]+"\t";

			for (int x=0;x<par;x++){  //
				line += IJ.d2s(res[y][x],4)+"\t";
			}
			if (tw==null) {
				title = "Surface Characterisation...";

				tw = new TextWindow(title,headings,line, width, (height));
			}
			else
			tw.append(line);
		}
	}

	void showInfo(String information){//, int n) {
			String infoTitle= "Information file";
			TextWindow info = new TextWindow(infoTitle, information, 400, 150);
	}
}