package active_contour;

import java.awt.*;
import java.awt.image.BufferedImage;
import java.io.*;
import java.awt.image.WritableRaster;
import java.awt.image.Raster;
import java.lang.Cloneable;
import java.lang.reflect.Array;
import java.util.ArrayList;
import java.util.Iterator;

import org.opencv.core.*;
import org.opencv.highgui.Highgui;
import org.opencv.imgproc.Imgproc;

import javax.imageio.ImageIO;
import javax.swing.ImageIcon;
import javax.swing.JFrame;
import javax.swing.JLabel;

/**
 * @version 25.4.2015
 * @author Branislav Hesko
 * This class calculates snake using active contour mdoel without edges
 */


public class ActiveContourWe {
	
	// Inicializacia premennych, niektore su uz finalne nastavene na urcitu hodnotu
	private ArrayList<Double> curve;
	private static final double lambda1 = -1;
	private static final double lambda2 = 1;
	private static final int sizeStructureElement = 5;
	private double stepT;
	private double stepH;
	private ImageProc img;
	private int numberOfPoints;
	private int numberOfIter;
	private double c1,c2;
	private double perimeter;
	private Mat subImg;
	private Mat image;
	private int matrixSizeY, matrixSizeX;
	private Mat phi;
	private int rowStart,rowEnd;
	private double mu;
	private double vu;
	private Mat result;
	private Mat result2;
	private int middleX;
	private int middleY;
	
	/*
	 * Konstruktor, ktory priradi lokalnym premennym predavane parametre z triedy main a tento konstruktor obsahuje cely vypocet
	 * metody.  
	 */
	public ActiveContourWe(ImageProc img,int middleX, int middleY, int numberOfIter, double perimeter, double stepH, double stepT, double mu, double vu , int rowStart, int rowEnd, int matrixSizeX, int matrixSizeY) {
		// Priradenie predavanych premennych a transoformacia obrazu na typ double
		this.img = img;
		this.image = new Mat(img.getNumRows(),img.getNumCols(), CvType.CV_64F);
		img.getMat().convertTo(image, CvType.CV_64F);
		this.numberOfIter = numberOfIter;
		this.perimeter = perimeter;
		this.matrixSizeX = matrixSizeX;
		this.matrixSizeY = matrixSizeY;
		this.rowStart = rowStart;
		this.rowEnd = rowEnd;
		this.stepH = stepH;
		this.stepT = stepT;
		this.mu = mu;
		this.vu = vu;
		this.middleX = middleX;
		this.middleY = middleY;
		
		// Inicializacia nulovych matic pre jednotlive premenne, phi - matica parametrickeho priestoru
		this.phi = new Mat(matrixSizeX,matrixSizeY,CvType.CV_64F);
		// Matica vyslednej kontury
		this.result = new Mat(matrixSizeX,matrixSizeY,CvType.CV_64F);
		// Pomocna matica
		this.result2 = new Mat(matrixSizeX,matrixSizeY,CvType.CV_64F);
		// Matica s originalnym vyrezanym obrazom
		this.subImg = new Mat(matrixSizeX,matrixSizeY,CvType.CV_64F);
		// Vypocet vyerzanej matice
		getSubMatrix().copyTo(subImg);
		// Privolanie inicializacie krivky
		inicialiseCurve();
		// Privolanie vyvoja krivky
		evolveCurve();
		// Vytvorenie binarnej masky vnutra krivky
		makeBinaryMask();
		
		// Ulozenie binarnej masky do suboru
		Highgui.imwrite("orig.png", result);
		// Dilatacia binarnej masky
		dilateImage();
		// Pomocna operacia
		result2.copyTo(result);
		// Erozia obrazu
		erodeImage();
		// Pomocne operacie a ulozenie do suboru
		result2.copyTo(result);
		Highgui.imwrite("dilero.png", result);
		
		// Detekcia KONTURY!
		findContour();
		result2.copyTo(result);
		// Ulozenie do suboru kontury
		Highgui.imwrite("hra.png", result);
		
	}
	/* 
	 * Inicializacia krivky
	 */
	public void inicialiseCurve() {
		int i,j;
		for(i = 0;i < matrixSizeX;i++) {
			for(j = 0;j <  matrixSizeY;j++) {
				// Naplnenie parametrickeho priestoru
				phi.put(i,j, (- Math.sqrt( Math.pow(j-middleX,2) + Math.pow(i-middleY,2)) + perimeter));
			}
		}
	}
	
	/*
	 * HLAVNA metoda - zabezpecuje vyvoj krivky
	 */
	public void evolveCurve() {
		int i,j,k,l,m;
		double term1,term2,term3;
		double result;
		double difC= 0;
		Mat phiZaloha = new Mat(matrixSizeX,matrixSizeY,CvType.CV_64F);
		int AD;
		// Zaciatok iterativneho procesu
		for(m = 0;m<numberOfIter;m++) {
			// Vypocet konstant c1 a c2 ako priemer hodnot mimo a vnutri krivky
			calculateMean();
			System.out.println("Hodnota konstanty c1: "+c1 + " Hodnota konstanty c2: " +c2);
			
			// Osetrenie skokovej divergencie, ak sa zmensi rozdiel konstant, ukonci sa proces a povodna hodnota phi sa
			// skopiruje do vysledku
			if(Math.abs(c1-c2) < difC) {
				phiZaloha.copyTo(phi);
				System.out.println("Iterativny vyvoj bol zastaveny, proces zacal divergovat");
				break;
			}else {
				// 0.5, pretoze nepatrne kolisanie rozdielu je povolene
				difC = Math.abs(c1-c2)-1;
				phi.copyTo(phiZaloha);
			}
			for(i = 1;i<matrixSizeX-1;i++) {
				for(j=1;j<matrixSizeY-1;j++) {
					// Iterativny vypocet nasledujuceho kroku, najlepsia volba je hodnota v ktorku m
					result = phi.get(i, j)[0];
					// Iterativne sa pocitaju vsetky casti rovnice z diplomovej prace 
					for(k = 0;k<3;k++) {
					term1 = (getDeltaPlusX(i,j)*result) / Math.sqrt(Math.pow((phi.get(i,j+1)[0] - phi.get(i,j-1)[0]),2)/4.0/stepH/stepH 
							+ Math.pow((getDeltaPlusX(i,j)*phi.get(i,j)[0]),2)/stepH/stepH )*mu/stepH/stepH*getDeltaMinusX(i,j);
					//System.out.println("Term1: " + Math.pow((phi.get(i,j+1)[0] - phi.get(i,j-1)[0]),2));
					term2 = (getDeltaPlusY(i,j)*result) / Math.sqrt(Math.pow((phi.get(i+1,j)[0] - phi.get(i-1,j)[0]),2)/4.0/stepH/stepH 
							+ Math.pow((getDeltaPlusY(i,j)*phi.get(i,j)[0]),2)/stepH/stepH )*mu/stepH/stepH*getDeltaMinusY(i,j);
					//System.out.println("Term2: " + Math.pow((getDeltaPlusX(i,j)*phi.get(i,j)[0]),2));

					term3 = -vu + lambda1*Math.pow((subImg.get(i, j)[0] - c1),2) + lambda2*Math.pow((subImg.get(i, j)[0] - c2),2);
	
					
					//System.out.println("Term3:  " + term2);
					// Vypocet hodnoty vysledku v jeho iterativnom vypocte s pouzitim heavisiideovej funkcie
					result = calculateHeaviside(i,j) * (term1+term2+(term3))*stepT + phi.get(i, j)[0];
					//result = calculateHeaviside(i,j) * (term3)*stepT + phi.get(i, j)[0];
					
					}
					//System.out.println("Result: " + result);
					// Uozenie do parametrickeho priestoru vyslednej hodnoty
					phi.put(i, j, result);
					
				}
			}
		}
	}
	
	/*
	 * Metoda vypoctu priemeru, c1 je priemer vnutri krivky, kde parameter je vacsi ako nula a c2 priemer mimo krivky,
	 * kde parameter mensi ako nula
	 */
	public void calculateMean() {
		double numberIn = 0,numberOut = 0;
		int i,j;
		double[] var,var1;

		c1 = 0;
		c2 = 0;
		for(i = 0; i<matrixSizeX;i++) {
			for(j = 0; j<matrixSizeY;j++) {
		
				var = phi.get(i,j);
				 if (var[0] > 0) {
					numberIn++;
					var1 = subImg.get(i,j);
					c1 = c1 + var1[0]; 
				 }
				 else if(var[0] < 0) {
					numberOut++;
					var1 = subImg.get(i,j);
					c2 = c2 + var1[0];  
				 }
			}
		}
		c1 = c1/numberIn;
		c2 = c2/numberOut;
	}
	
	/* Ziskanie vyrezaneho obrazu na spracovanie
	 * 
	 */
	public Mat getSubMatrix() {
		return image.submat(rowStart, rowStart + matrixSizeX,rowEnd , rowEnd + matrixSizeY);
	}
	/*
	* Nevyuzita metoda pre vypocet dirakovej funkcie
	*/
	public double calculateDirac(int i, int j) {
		double[] var;
		double delta;
		var = phi.get(i, j);
 		delta = 1.0/(var[0]* var[0] + 1);
 		return delta;
	}

	/*
	 * Vypocet alternativne zadanej heavisideovej funkcie
	 */
	public double calculateHeaviside(int i,int j) {
		double[] var;
		double k = 1;
		double delta;
		var = phi.get(i, j);
		delta = 1.0/2.0*(1.0 + 2.0/Math.PI*Math.atan(var[0])/k);
		//System.out.println("delta  " +delta);
		return delta;
	}
	
	/*
	 * Zisk doprednej diferencie v osi x
	 */
	public double getDeltaPlusX(int i, int j) {
			double[] var, var1;
			var = phi.get(i+1,j); 
			var1 = phi.get(i,j);
			//System.out.println(var[0] + "   "  +var1[0] );
		return  (var[0] - var1[0]);
	}
	
	/*
	 * Zisk doprednej diferencie v osi y
	 */
	public double getDeltaPlusY(int i, int j) {
		double[] var, var1;
		var = phi.get(i,j+1); 
		var1 = phi.get(i,j);
	return  (var[0] - var1[0]);
	}
	
	/*
	 * Zisk spatnej diferencie v osi x
	 */
	public double getDeltaMinusX(int i, int j) {
		
		double[] var, var1;
		var = phi.get(i,j); 
		var1 = phi.get(i-1,j);
	return  (var[0] - var1[0]);
	}
	
	/*
	 * Zisk spatnej diferencie v osi y
	 */
	public double getDeltaMinusY(int i, int j) {
		
		double[] var, var1;
		var = phi.get(i,j); 
		var1 = phi.get(i-1,j);
	return  (var[0] - var1[0]);
	}
	
	/*
	 * Zisk cisla riadka pri pouziti vektoroveho vyjadrenia matice
	 */
	public int getXaxe(Object o) {
		int var = ((int) o) / img.getNumCols();	
		return var;
	}
	
	/*
	 * Zisk cisla stlpca pri pouziti vektoroveho zapisu matice
	 */
	public int getYaxe(Object o) {
		int var = ((int) o) / img.getNumCols();
		var = ((int) o)  - (var*img.getNumCols());
		return var;
	}
	
	/*
	 * Ziska binarnej masky ako vnutra obrazu, kedy je parameter kladny
	 */
	public void makeBinaryMask() {
		int i,j;
		for(i = 0;i<matrixSizeX;i++) {
			for(j = 0;j<matrixSizeY;j++) {
				if(phi.get(i, j)[0] > 0) {
					result.put(i, j, 255);
				} else {
					result.put(i, j, 0);
				}
			}
		}
	}
	
	/*
	 * Dilatacia obrazu
	 */
	public void dilateImage() {
		int i,j,k,l;
		int var;
		int var2 = sizeStructureElement/2;
		System.out.println("var2: " +var2);
		for(i= var2;i < matrixSizeX-var2;i++) {
			for(j = var2; j< matrixSizeY-var2; j++) {
				var = 0;
				for(k = 0;k < sizeStructureElement; k++) {
					for(l = 0;l<sizeStructureElement; l++) {
						if(result.get((i-var2 + k), (j- var2 +l))[0] > 0) {
							var++;
						}
					}
				}
				if(var > 0) {
					result2.put(i,j,255);
				} else {
					result2.put(i,j,0);
				}
				
			}
		}
	}
	
	/*
	 * Erozia obrazu
	 */
	public void erodeImage() {
		int i,j,k,l;
		int var;
		int sizeStructureElement2 = sizeStructureElement + 2;
		int var2 = sizeStructureElement2/2;
		
		for(i= var2;i < matrixSizeX-var2;i++) {
			for(j = var2; j< matrixSizeY-var2; j++) {
				var = 0;
				for(k = 0;k < sizeStructureElement2; k++) {
					for(l = 0;l<sizeStructureElement2; l++) {
						if(result.get((i-var2 + k), (j- var2 +l))[0] > 0) {
							var++;
						}
					}
				}
				if(var == (sizeStructureElement2*sizeStructureElement2)) {
					result2.put(i,j,255);
				} else {
					result2.put(i,j,0);
				}
				
			}
		}
	}
	
	/*
	 * Detekcia kontury
	 */
	public void findContour() {
		int i,j,k,l;
		int var;
		int var2 = 3/2;
		
		
		for(i= var2;i < matrixSizeX-var2;i++) {
			for(j = var2; j< matrixSizeY-var2; j++) {
				var = 0;
				for(k = 0;k < 3; k++) {
					for(l = 0;l<3; l++) {
						if(result.get((i-var2 + k), (j- var2 +l))[0] > 0) {
							var++;
						}
					}
				}
				if(var == (3*3) || var < 5) {
					result2.put(i,j,0);
				} else {
					result2.put(i,j,255);
				}
				
			}
		}
	}
	/*
	 * Vrati vyslednu konturu
	 */
	public Mat getContour() {
		return result;
	}
	
	/*
	 * Vrati maticu parametrickeho priestoru
	 */
	public Mat getParametric() {
		return phi;
	}
	
	
}

