import ij.*;
import ij.plugin.*;
import ij.process.*;
import ij.gui.*;
import java.awt.*;
import java.util.Vector;
import java.awt.event.*;
import java.io.*;
import java.lang.Math;
import ij.plugin.filter.PlugInFilter;

public class CoLocalization_ implements PlugInFilter {

	private static final int MAX = 2000;
	private int num1 = 0;
	private int num2 = 0;
	private final Vector listPoints1 = new Vector(0,16);
	private final Vector listClus1 = new Vector(0,16);
	private int clust_n1=0;
	private final Vector listPoints2 = new Vector(0,16);
	private final Vector listClus2 = new Vector(0,16);
	private int clust_n2=0;
	private double thresh_d1 = 105;
	private double thresh_d2 = 50;
	private double thresh_d3 = 52;
	private int co1=0;
	private int co2=0;
	private int max1=0;
	private int max2=0;
	private final Vector listClusPair = new Vector(0,16);

	public int setup(String ar, ImagePlus imp){
		return DOES_8G;
	}

	public void run(ImageProcessor ip){
		GenericDialog dia1 = new GenericDialog("Set Threshold", IJ.getInstance());
		dia1.addNumericField("Threshold (10nm)", thresh_d1, 0);
		dia1.addNumericField("Threshold (5nm)",thresh_d2,0);
		dia1.addNumericField("Threshold (co)",thresh_d3,0);
		dia1.showDialog();

		if(dia1.wasCanceled()) return;
		if(dia1.invalidNumber()) {
				IJ.showMessage("Error", "Invalid input Number");
				return;
				}

		thresh_d1 = dia1.getNextNumber();
		thresh_d2 = dia1.getNextNumber();
		thresh_d3 = dia1.getNextNumber();
		IJ.showMessage("Please open 10nm file");
		openFile(listPoints1);
		IJ.showMessage("Please open 5nm file");
		openFile(listPoints2);
		IJ.write("**Threshold for 10nm: "+thresh_d1+"\n");
		IJ.write("**Threshold for 5nm: "+thresh_d2+"\n");
		IJ.write("**Threshold for CoLocalization: "+thresh_d3+"\n");
		IJ.write("**Number of 10 nm particles: "+num1+"\n**Number of 5nm particles: "+num2+"\n");
		calculate_neighbors(listPoints1);
		calculate_neighbors(listPoints2);
		cal_cluster(listPoints1,listClus1);
		cal_cluster(listPoints2,listClus2);
		IJ.write("**Number of clusters of 10nm particles: "+clust_n1+"\n**Number of clusters of 5nm particles: "+clust_n2+"\n");

		find_co_cluster();

        ImageProcessor ip2=ip.duplicate();

		ParticleCluster pc;
		ip2.setLineWidth(2);

		for(int i=0;i<clust_n1;i++){
			pc=(ParticleCluster)listClus1.elementAt(i);
			ip2.drawLine(pc.min_x-20,pc.min_y-20,pc.max_x+20,pc.min_y-20);
			ip2.drawLine(pc.min_x-20,pc.min_y-20,pc.min_x-20,pc.max_y+20);
			ip2.drawLine(pc.min_x-20,pc.max_y+20,pc.max_x+20,pc.max_y+20);
			ip2.drawLine(pc.max_x+20,pc.max_y+20,pc.max_x+20,pc.min_y-20);
		}

		for(int j=0;j<clust_n2;j++){
			pc=(ParticleCluster) listClus2.elementAt(j);
			int r,x,y;
			x=(pc.min_x+pc.max_x)/2;
			y=(pc.min_y+pc.max_y)/2;
			int r2=(pc.max_x-pc.min_x)*(pc.max_x-pc.min_x)+(pc.max_y-pc.min_y)*(pc.max_y-pc.min_y);
			if(r2==0) r=20;
			else r=((int) Math.pow(r2,0.5))/2 +10;
			drawCircle(r,x,y,ip2);
		}

        new ImagePlus("Clutering Results",ip2).show();

		IJ.write("**co1 = "+co1+" co2 = "+co2+"\n");
		IJ.write("**10 nm particles : "+(double)co1/num1+" 5nm particles : "+(double)co2/num2+"\n");
		IJ.write(max1+"\n"+max2+"\n");
		int k;
		int [] array1 = new int[max1+1];
		int [] array2 = new int[max2+1];
		for(k=0;k<max1+1;k++) {
			array1[k]=0;
		}
		for(k=0;k<max2+1;k++) {
			array2[k]=0;
		}


		//IJ.write("**10nm particle clusters:");
		for(k=0;k<listClus1.size();k++){
			pc= (ParticleCluster) listClus1.elementAt(k);
			if(!pc.colocal) array1[pc.getNum()]++;
		}
		int nn10=0;
		for(k=1;k<max1+1;k++) {
			if(!(array1[k]==0)) nn10++;
		}
		IJ.write(nn10+"\n");
		for(k=1;k<max1+1;k++) {
			if(!(array1[k]==0)) IJ.write(k+" "+array1[k]+"\n");
		}

		//IJ.write("**5nm particle clusters:");
		for(k=0;k<listClus2.size();k++){
			pc= (ParticleCluster) listClus2.elementAt(k);
			if(!pc.colocal) array2[pc.getNum()]++;
		}
		int nn5=0;
		for(k=1;k<max2+1;k++) {
			if(!(array2[k]==0)) nn5++;
		}
		IJ.write(nn5+"\n");
		for(k=1;k<max2+1;k++) {
			if(!(array2[k]==0)) IJ.write(k+" "+array2[k]+"\n");
		}

		//IJ.write("**coclustering:");
		combineClusterPair();

/*
		for(k=0;k<listClusPair.size();k++){
			ParticleClusterPair ppp= (ParticleClusterPair) listClusPair.elementAt(k);
			IJ.write(ppp.pc1.getID()+" "+ppp.pc2.getID()+" "+ppp.duplicatePC1);
		}
*/


	}

	void drawCircle(int r, int x, int y, ImageProcessor ip){
		double radius = (double) r;
		for(int var=1;var<=2;var++){
			for (double counter = 0; counter < 10; counter = counter + 0.001) {
				double xx = Math.sin(counter) * (radius+var) + x;
				double yy = Math.cos(counter) * (radius+var) + y;
				ip.drawPixel((int)xx, (int)yy);
			}
		}
	}

	void find_co_cluster(){
		ParticleCluster pc1, pc2;
		for(int i=0;i<clust_n1;i++){
			pc1 = (ParticleCluster) listClus1.elementAt(i);
			for(int j=0;j<clust_n2&&(!pc1.colocal);j++){

				pc2 = (ParticleCluster) listClus2.elementAt(j);
				if(pc1.IsCoLocalization(pc2,thresh_d3)) {
					co1+=pc1.getNum();
					//IJ.write(pc1.getID()+" "+pc2.getID());
					//IJ.write("find cocluster of "+pc1.getNum()+" 10 nm particles and "+pc2.getNum()+" 5nm particles\n");
					pc1.colocal=true;
					pc2.checked=true;
					ParticleClusterPair pcp = new ParticleClusterPair();
					pcp.pc1=pc1;
					pcp.pc2=pc2;
					listClusPair.addElement(pcp);
				}
			}
		}
		//IJ.write("++++++++++++++++++++++++++++++++");
		for(int i=0;i<clust_n2;i++){
			pc1 = (ParticleCluster) listClus2.elementAt(i);
			for(int j=0;j<clust_n1&&(!pc1.colocal);j++){

				pc2 = (ParticleCluster) listClus1.elementAt(j);
				if(pc1.IsCoLocalization(pc2,thresh_d3)) {
					co2+=pc1.getNum();
					if(!pc1.checked) {
						//IJ.write(pc2.getID()+" "+pc1.getID());
						//IJ.write("find cocluster of "+pc2.getNum()+" 10 nm particles and "+pc1.getNum()+" 5nm particles\n");
						ParticleClusterPair pcp = new ParticleClusterPair();
						pcp.pc1=pc2;
						pcp.pc2=pc1;
						listClusPair.addElement(pcp);
					}
					pc1.colocal=true;

				}
			}
		}
	}

	void combineClusterPair(){
		ParticleClusterPair pp, pp1;
		int numPC2;

		int max_5co=2; //max number of 5nm clusters with the same 10nm cluster
		int temp;
		for(int i=0;i<listClusPair.size();i++){
			pp = (ParticleClusterPair) listClusPair.elementAt(i);
			temp=1;
			for(int j=0;j<listClusPair.size();j++){
				pp1 = (ParticleClusterPair) listClusPair.elementAt(j);
				if(pp.pc1.getID()==pp1.pc1.getID()){
					temp++;
				}
			}
			//IJ.write(temp+" ");
			if(temp>max_5co) max_5co=temp;
		}
		int [] [] array3 = new int[max1+1][max2*max_5co];
		for(int i=0;i<max1+1;i++){
			for(int j=0;j<max2*max_5co;j++){
				array3[i][j]=0;
			}
		}

		for(int i=0;i<listClusPair.size();i++){
			pp = (ParticleClusterPair) listClusPair.elementAt(i);
			numPC2=0;
			if(!pp.duplicatePC1){
				for(int j=i+1;j<listClusPair.size();j++){
					pp1 = (ParticleClusterPair) listClusPair.elementAt(j);
					if((!pp1.duplicatePC1)&&(pp.pc1.getID()==pp1.pc1.getID())){
						numPC2+=pp1.pc2.getNum();
						pp1.duplicatePC1=true;
					}
				}
			}
			if(!pp.duplicatePC1) {
				numPC2+=pp.pc2.getNum();
				if(numPC2>=max2*max_5co) IJ.write("** ERROR** 5nm Overflow\n");
				array3[pp.pc1.getNum()][numPC2]++;
			}
		}

		int nnc=0;
		for(int i=0;i<max1+1;i++){
			for(int j=0;j<max2*max_5co;j++){
				if(!(array3[i][j]==0)) nnc++;
			}
		}
		IJ.write(nnc+"\n");
		int xmm=0;
		int ymm=0;
		for(int i=0;i<max1+1;i++){
			for(int j=0;j<max2*max_5co;j++){
				if(!(array3[i][j]==0)) {
					IJ.write(i+" "+j+" "+array3[i][j]+"\n");
					if(i>xmm) xmm=i;
					if(j>ymm) ymm=j;
				}
			}
		}
		IJ.write(xmm+"\n");
		IJ.write(ymm+"\n");
	}




	void openFile(Vector v){
		Frame f = new Frame();
		FileDialog fd = new FileDialog(f,"point list",FileDialog.LOAD);
		fd.setVisible(true);
		String path = fd.getDirectory();
		String filename = fd.getFile();
		if((path==null)||(filename==null)){return;}
		try{
			FileReader fr = new FileReader(path+filename);
			BufferedReader br = new BufferedReader(fr);
			String line;
			String pString;
			String xString;
			String yString;
			int separatorIndex;
			int x;
			int y;
			if ((line = br.readLine()) == null) {
				fr.close();
				return;
			}
			while ((line = br.readLine()) != null) {
				line = line.trim();
				//if(line.length()==0) return; //empty line
				separatorIndex = line.indexOf(' ');
				if (separatorIndex == -1) {
					fr.close();
					IJ.error("Invalid file");
					return;
				}
				line = line.substring(separatorIndex);
				line = line.trim();
				separatorIndex = line.indexOf(' ');
				if (separatorIndex == -1) {
					fr.close();
					IJ.error("Invalid file");
					return;
				}
				xString = line.substring(0, separatorIndex);
				xString = xString.trim();
				line = line.substring(separatorIndex);
				line = line.trim();
				separatorIndex = line.indexOf(' ');
				if (separatorIndex == -1) {
					separatorIndex = line.length();
				}
				yString = line.substring(0, separatorIndex);
				yString = yString.trim();
				//IJ.write(xString+" "+yString);
				x = Integer.parseInt(xString);
				y = Integer.parseInt(yString);
				addPoint(x, y, v);
			}
			fr.close();
		}catch(FileNotFoundException e){
			IJ.error("file not found exception");
		}catch(IOException e){
			IJ.error("IO exception");
		}
	}

	public void addPoint(int x, int y, Vector v){
		int num=0;
		if(v==listPoints1) num=num1;
		if(v==listPoints2) num=num2;
		if(num<MAX){
			Particle p = new Particle(x,y,0);
			v.addElement(p);
			if(v==listPoints1) num1++;
			if(v==listPoints2) num2++;
		}else{IJ.error("Maxium number of points reached");}
	}

	public void calculate_neighbors(Vector v){

		double d,d1;
		Particle p1;
		Particle p2;
		double thresh_d=105;
		if(v==listPoints1) thresh_d=thresh_d1;
		if(v==listPoints2) thresh_d=thresh_d2;

		for(int k=0;k<v.size();k++){
			p1=(Particle) v.elementAt(k);

			for(int j=0;j<v.size();j++){
				if(j!=k){
					p2=(Particle) v.elementAt(j);
					d1 = (p2.x-p1.x)*(p2.x-p1.x) + (p2.y-p1.y)*(p2.y-p1.y);
					d = Math.pow((double)d1,0.5);
					if(d<thresh_d) p1.addNeighbor(p2);
				}
			}
		}
	}


	public void cal_cluster(Vector v1, Vector v2){
		Particle p1,p2;
		int clust_id=1;
		ParticleCluster pc;
		double dd,d;
		int limit;


		for(int i=0;i<v1.size();i++){
			p1=(Particle) v1.elementAt(i);

			if(p1.cnum==0){
			pc= new ParticleCluster(clust_id);
			pc.addParticle(p1);
			p1.cnum=clust_id;
			pc.findAllParticles();


            v2.addElement(pc);
            if(v2==listClus1) {
				clust_n1++;
				if(pc.getNum()>max1) max1=pc.getNum();
			}
            if(v2==listClus2) {
				clust_n2++;
				if(pc.getNum()>max2) max2=pc.getNum();
			}
			clust_id++;
		}



		}
	}


}

class Particle {
	public int x;
	public int y;
	public int cnum;
	public int nn;
	public Vector neighbors = new Vector(0,16);

	Particle(int a,int b, int c){
		x=a;
		y=b;
		cnum=c;
		nn=0;
	}

	public void addNeighbor(Particle p){
		neighbors.addElement(p);
		nn++;
	}

	public void print(){
		IJ.write("X= "+x+", Y= "+y+"Cluster #:"+cnum+"\n");
	}
}


class ParticleCluster {
	private Vector particles = new Vector(0,16);
	private Vector PCNeighbors = new Vector(0,16);
	private int num_p =0;
	private int neigb= 0;
	public int max_x = 0;
	public int max_y = 0;
	public int min_x = 0;
	public int min_y = 0;
	private int cluster_id = 0;
	public boolean colocal = false;
	public boolean checked = false;

	public void addParticle(Particle p){

		if(num_p==0) {
			max_x=p.x;
			min_x=p.x;
			max_y=p.y;
			min_y=p.y;
		}

		p.cnum=cluster_id;
		particles.addElement(p);
		addNeighbor(p);
		num_p++;
		if(p.x>max_x) max_x=p.x;
		if(p.y>max_y) max_y=p.y;
		if(p.x<min_x) min_x=p.x;
		if(p.y<min_y) min_y=p.y;
	}

	public double distance(Particle p,int x, int y){
		double d1;
		d1 = (p.x-x)*(p.x-x) + (p.y-y)*(p.y-y);
		return Math.pow((double)d1,0.5);
	}

	public double getMinD(int x,int y){
		if(num_p==0) return -1.0;

		double result,temp;
		Particle p1,p;

		p1=(Particle) particles.elementAt(0);
		result = distance(p1,x,y);

		for(int i=1;i<num_p;i++){
			p=(Particle) particles.elementAt(i);
			temp = distance(p,x,y);
			if(temp<result) result=temp;
		}

		return result;
	}



	public void setID(int i){
		cluster_id=i;
	}

	public int getID(){
		return cluster_id;
	}

	public int getNum(){
		return num_p;
	}

	public void print(){
		IJ.write("Cluster #:"+cluster_id+" Number of particles :"+num_p+"\n");
		IJ.write(" x from "+min_x+" to "+max_x+" y from "+min_y+" to "+max_y+"\n");
	}

	ParticleCluster (int id){
		cluster_id=id;
	}


	private void addNeighbor(Particle p){
		Particle p1;
		for(int i=0;i<p.nn;i++){
			p1=(Particle) p.neighbors.elementAt(i);
			if(!InCluster(p1)) {
				PCNeighbors.addElement(p1);
				neigb++;
			}
		}
	}

	private int CloseNeighbor(){
		int result;
		double d1,d2;
		Particle p1,p2;

		p1=(Particle) PCNeighbors.elementAt(0);
		d1=getMinD(p1.x,p1.y);
		result=0;
		for(int i=1;i<neigb;i++){
			p2=(Particle) PCNeighbors.elementAt(i);
			d2=getMinD(p2.x,p2.y);
			if(d2<d1) result=i;
		}

		return result;
	}


    public void findAllParticles(){
		int next;
		int limit=0;
		Particle p;
		while(neigb>0&&limit<1000){
			next=CloseNeighbor();
			p=(Particle) PCNeighbors.elementAt(next);
			PCNeighbors.removeElementAt(next);
			neigb--;
			if(!InCluster(p)) addParticle(p);
			limit++;

		}
	}

	public boolean InCluster(Particle p){
		boolean result=false;
		Particle p1;
		for(int i=0;i<num_p;i++){
			p1=(Particle) particles.elementAt(i);
			if(p1==p) result=true;
		}

		return result;
	}

	public boolean IsCoLocalization(ParticleCluster pc, double d){
		boolean result=false;
		Particle p1, p2;
		for(int i=0;i<num_p;i++){
			for(int j=0;j<pc.num_p;j++){
				p1=(Particle) particles.elementAt(i);
				p2=(Particle) pc.particles.elementAt(j);
				double dd = distance(p1,p2.x,p2.y);
				if(dd<d) {
					result = true;
					break;
				}
			}
		}

		return result;
	}



}

class ParticleClusterPair {
	public ParticleCluster pc1;
	public ParticleCluster pc2;
	boolean duplicatePC1 = false;
}