/**** Arvonn Tully 3/19/2004
      This code is public domain, but I would appreciate if this code is used
      you would acknowlegde my contributions, and the contributions of 
      Wayne Rasband for his help, and example code. 

      Purpose:
      Plugin for ImageJ to calculate a correlation coefficent between images

      Requirements:
      Input: 
      1) One stack of 16bpp Grayscale Image, minimum 2 images, no maximum, with ROI 

      Output:
      1) a list of (n-1) doubles, 0>=X>=1(in theory), where n = number of images, where each number is the 
      pixel by pixel correlation coefficent of the user selected ROI between adjacent images
      2) a list of (n-1) vectors, where n = number of images, where each vector is 4 numbers,
      the average correlation coeffficent of the user selected ROI between images, number of images
      used to make the average, Lag, and Standard Error for each set. 

      Other Goals:
      1) Code reusablity for development of other ImageJ plugins - especially for image registration.
      2) Intuitive user interface, with helpful error messages. 
      3) Accuracy, ie, same answers as Innovision's Isee Software, which this plugin is intended as 
      a free and portable replacement for the Correlation Coefficent output from the Nanotracker pnode. 
*/

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


public class CC_calc implements PlugIn {
    private static int index, lag, start;
    private static boolean autocorrelate, show_stats, show_plot;
    private static ImageStack my_stack;
    private static Roi my_roi; 
    private static Rectangle r;
    private static ImageProcessor mask;
    
    
    public void run(String arg){
	if (IJ.versionLessThan("1.32c")){
	    IJ.showMessage("Please upgrade ImageJ to at least version 1.32c");
	    return;
	}
	if (input()) 
	    calc();
	else IJ.showMessage("Please Provide Valid Input");
    }//run 

    public boolean input(){
	/*
	  Input:
	  1. Get Valid Imagestack(16bit grayscale) 
	  2. Get Valid ROI (!= null)
	  3. Correlation(false) or Autocorrelation(true)
	 
	  Output:
	  Returns true on sucessful setting of following private variables
	  1. ImageStack my_stack
	  2. boolean autocorrelate
	  3. Roi my_roi
	*/
   
        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 (index>=titles.length)index = 0;

        GenericDialog gd = new GenericDialog("Correlation Coefficent Calculator");
        gd.addChoice("Image Stack :", titles, titles[index]);
	gd.addNumericField("Start :", 1, 0);
	gd.addNumericField("Lag :", 1, 0);

	gd.addCheckbox("Autocorrelate : ", autocorrelate);
	gd.addCheckbox("Show Autocorrelation Statistics: ", show_stats);
	gd.addCheckbox("Show Stats plot: ", show_plot);
        gd.showDialog();

        if (gd.wasCanceled()) {
	    return false;
	}
	
        index = gd.getNextChoiceIndex();
        ImagePlus sliceimg = WindowManager.getImage(wList[index]);

        if (sliceimg.getType()!=sliceimg.GRAY16) {
            IJ.showMessage("Image Correlator", "The stack must be 16-bit grayscale.");
            return false;
        }

	//set values
	start = (int)gd.getNextNumber();
	lag = (int)gd.getNextNumber();
	autocorrelate = gd.getNextBoolean();
	show_stats = gd.getNextBoolean();
	show_plot = gd.getNextBoolean();
	my_roi = sliceimg.getRoi();
        my_stack = sliceimg.getStack();
	r=my_roi.getBoundingRect(); 
	mask=my_roi.getMask();

        return true;
    }//input


    public void calc(){
        if(start>my_stack.getSize()){
	    IJ.showMessage("Please pick a valid starting point.");
	    return;
	}
	if (autocorrelate){
	    /*
	      Basic idea of autocorrelation:
	      obtain average CC of a certain lag for all images, average that lag.
	      repeat for all lags from 1 to size of array. 
	      return average, std. error, and n for each lag. 

	    */
	    
	    //setup values for w or w/o stats
	    int n, auto_lag; 

	    n = 0; auto_lag =0;
	    ResultsTable rt=new ResultsTable();
 
	    if (show_stats){
		//these values only needed for showing statistics
		Vector current_doubles = new Vector();
		double avg, std_err;
		avg = 0; std_err =0;

		rt.setHeading(1,"Average");
		rt.setHeading(2,"Standard Error");
		rt.setHeading(3,"N");
		rt.setHeading(4,"Lag");	    	    

		rt.disableRowLabels();
		rt.setPrecision(6);

		/*
		//Vector Notes 
		void current_doubles.addElement(obj);//add obj, increase size by 1
		void current_doubles.setSize(int);//increase size eg. allocated memory
		void current_doubles.trimToSize();//reduces memory allocated
		void current_doubles.removeAllElements();//emptys the vector
		void current_doubles.copyInto(Object [] anarray);//copys contents into an array.
		*/
		n=my_stack.getSize();
		double [] avgs = new double[n];
		float [] errs = new float[n];
		double [] xvals = new double[n];
		for(auto_lag=1; auto_lag<=(n-1);auto_lag++){
		    current_doubles.removeAllElements();
		    //n=0;
		    for(int i=1; i<=(my_stack.getSize()-auto_lag);i++){
			//n++;
			current_doubles.addElement(new Double(correlation_coefficent(my_stack.getProcessor(i),my_stack.getProcessor(i+auto_lag))));
		    }//rof i
		
		    std_err=standard_error(current_doubles);
		    avg=average(current_doubles);
		
		    //output to Results table. 
		    rt.incrementCounter();//have to do this first
		    rt.addValue(1,avg);
		    rt.addValue(2,std_err);
		    rt.addValue(3,(double)current_doubles.size());
		    rt.addValue(4,auto_lag);
		    
		    avgs[auto_lag]=avg;
		    errs[auto_lag]=(float)std_err;
		    xvals[auto_lag]=auto_lag;
		}//rof auto_lag
	    
		rt.show("Autocorrelation Results");
		if(show_plot){
		    String plot_title="Stats Plot";
		    String xLabel="Lag";
		    String yLabel="CC value";
		    PlotWindow stats_plot = new PlotWindow (plot_title, xLabel, yLabel, xvals, avgs);
		    stats_plot.addErrorBars(errs);
		    
		    double xMin = 1;
		    double xMax = n;
		    double yMin = avgs[n-1] - 0.01;
		    double yMax = 1.0;
		    stats_plot.setLimits(xMin, xMax, yMin, yMax);
		    stats_plot.draw();
		    
		    return;
		}
		return;
	    }
	    else{ // if (autocorrelate && !show_stats) - basically this out puts all the raw data from the CC
		rt.setPrecision(6);
		double currentvalue = 0;
		
		for(auto_lag=1; auto_lag<=(my_stack.getSize()-1);auto_lag++){
		    rt.incrementCounter();//increments the row
		    for(int i=1; i<=(my_stack.getSize()-auto_lag);i++){
			rt.setHeading(i, Integer.toString(i));
			currentvalue=correlation_coefficent(my_stack.getProcessor(i), my_stack.getProcessor(i+auto_lag));
			rt.addValue(i,currentvalue);
		    }//rof i
		}//rof auto_lag
	    
		rt.show("Autocorrelation Results");//shows table with title as label
		return;
	    }// not show_stats

	}//fi

	if (!autocorrelate){
	    /*
	    // this code snippet assumes lag == 1
	    for (int i = 1; i <= (my_stack.getSize()-1); i++){
	    IJ.write(Double.toString(correlation_coefficent(my_stack.getProcessor(i),my_stack.getProcessor(i+1))));
	    }
	    */

	    //lag==0 is special case, otherwise, lag is the distance between start and current
	    if(lag==0){
		for (int i = start; i <= my_stack.getSize(); i++){
		    IJ.write(Double.toString(correlation_coefficent(my_stack.getProcessor(start),my_stack.getProcessor(i))));
		}
	    }else {
		for (int i = start; i<=(my_stack.getSize()-lag); i++){
		    IJ.write(Double.toString(correlation_coefficent(my_stack.getProcessor(i),my_stack.getProcessor(i+lag))));
		}//rof
	    }//esle
	    return;
	}//fi
	
	return;
    }//calc



    public double standard_error(Vector my_set){
	double mean = average(my_set);//mean of set
	double var = 0; //variance
	double topsum=0;
	double curr = 0; // current 

	int n=my_set.size();
	if (n<3){
	    return 0; //std_err not valid for n < 3;
	}

	Double my_obj=new Double(0);

	for (int j=0;j<n;++j){
	    my_obj=(Double)my_set.get(j);
	    curr = my_obj.doubleValue() - mean;
	    topsum+=curr*curr;
	}
	var = topsum/(n-1);
	double answer= Math.sqrt(var);
	return(answer/Math.sqrt(n));
	//return(Math.sqrt(var));
    }




    /*
      //this one is wrong also
    public double standard_error(Vector my_set){
	double mean = average(my_set);//mean of set
	double var = 0; //variance
	double topsum=0;
	double curr = 0; // current 

	int n=my_set.size();

	Double my_obj=new Double(0);

	for (int j=0;j<n;++j){
	    my_obj=(Double)my_set.get(j);
	    curr = my_obj.doubleValue() - mean;
	    topsum+=curr*curr;
	}
	var = topsum/(n-1);
	//double answer= Math.sqrt(var);
	//return(answer);
	return(Math.sqrt(var));
    }
    */

    /*
      //this is closer than above, but still wrong
    public double standard_error(Vector my_set){
	double mean = average(my_set);
	double var = 0;
	double eps = 0;
	
	int n=my_set.size();
	Double my_obj=new Double(0);

	for (int j=0;j<n;++j){
	    my_obj=(Double)my_set.get(j);
	    double d = my_obj.doubleValue() - mean;
	    eps+=d;
	    var+=d*d;
	}
	var = (var - eps * eps/n)/(n-1);
	double stdev = Math.sqrt(var);
	double answer = (stdev/Math.sqrt(n));
	return(answer);
    }//std_error   
    */

    /*
      int size = vector.size();
      for( int index = 0; index < size; ++index )
      elementAt(index);
      }

    */

    public double average(Vector my_set){
	double answer =0;
	Double my_obj=new Double(0);
	int size =my_set.size();
	for(int i=0;i<size;++i){
	    my_obj=(Double)my_set.get(i);
	    answer+=my_obj.doubleValue();
	}
	return (answer/my_set.size());
    }//average
   
    public double correlation_coefficent(ImageProcessor ip1, ImageProcessor ip2){
	/*
	  Correlation Coefficent equals:
	  {N*Sum(XY)-Sum(X)*Sum(Y)}/SQRT({N*Sum(X*X)-Sum(X)*Sum(X)}*{N*Sum(Y*Y)-Sum(Y)*Sum(Y)})
	  Correlation Index equals:
	  1.0 - SQRT(SUM((Y-X)*(Y-X))/SUM(X*X))
	  
	  This function supports an  roi
	*/
	
	//if mask is not null, then it is a non-rectangular roi
	if(mask!=null){
	    double X,Y,N;
	    double CC,sumX, sumY, sumXX, sumXY, sumYY;
	    //init values
	    X=0;Y=0;N=0;CC=0;
	    sumX=0;sumY=0;sumXX=0;sumXY=0;sumYY=0;

	    for (int x=0;x<r.width;x++){
		for(int y=0;y<r.height;y++){
		    if (mask.getPixel(x,y)!=0){ //requires v 1.32c
			X=(double)ip1.getPixelValue(x+r.x,y+r.y);
			Y=(double)ip2.getPixelValue(x+r.x,y+r.y);
			N++;
			//for top
			sumXY+=X*Y; sumX+=X;sumY+=Y;
			//for bottom
			sumXX+=X*X;sumYY+=Y*Y;
		    }//mask
		}//rof height
	    }//rof width
	    double top = (N*sumXY)-(sumX*sumY);
	    double bottom = Math.sqrt(((N*sumXX)-(sumX*sumX))*((N*sumYY)-(sumY*sumY)));
	    CC = top/bottom;
	    
	    return CC;	
	}//fi my_roi exists
       
	//if mask is null, then it is a rectangular selection and we use the "r" as our roi. 
	if(mask==null&&r!=null){
	    double X,Y,N;
	    double CC,sumX, sumY, sumXX, sumXY, sumYY;
	    //init values
	    X=0;Y=0;N=0;CC=0;
	    sumX=0;sumY=0;sumXX=0;sumXY=0;sumYY=0;

	    for (int x=0;x<r.width;x++){
		for(int y=0;y<r.height;y++){
		    X=(double)ip1.getPixelValue(x+r.x,y+r.y);
		    Y=(double)ip2.getPixelValue(x+r.x,y+r.y);
		    N++;
		    //for top
		    sumXY+=X*Y; sumX+=X; sumY+=Y;
		    //for bottom
		    sumXX+=X*X;sumYY+=Y*Y;
		}//rof height
	    }//rof width
	    double top = (N*sumXY)-(sumX*sumY);
	    double bottom = Math.sqrt(((N*sumXX)-(sumX*sumX))*((N*sumYY)-(sumY*sumY)));
	    CC = top/bottom;
	    return CC;
	}//fi
	
	return 0;
    }//run_calc
}//cc_calc class