/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */
package kubovir;

import flanagan.analysis.Regression;
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.FileReader;
import java.io.IOException;
import java.io.OutputStreamWriter;
import java.io.Writer;
import java.util.ArrayList;
import java.util.List;

/**
 *
 * @copyright Monika Lucia Jakubcová
 */
public class KuboVir { 
    private final int leftBoundary = 100; //left boundary of cylinder tank
    private final int rightBoundary = 1030; //right boundary of cylinder tank
    private int middleOfContainer = leftBoundary + (rightBoundary - leftBoundary) / 2 ; //middle of cylinder tank
    private int middleOfWhirl; //middle of vortex
    
    private final List<List<Pixel>> intensities = new ArrayList<>();
    private List<Pixel> transitions = new ArrayList<>();
    
    private double yscale; //LORENTZIAN A0
    private double mean;  //LORENTZIAN mu
    private double gamma; //LORENTZIAN gamma
    private double sigma; //if GAUSSIAN
    
    private final int numberOfImages = 70; //number of images in the folder to process
    private int surfaceY; //water surface
    private int bottomY; //vortex peak
    
    private double searchLimit; //limit height for searching
    private double intensityTresholdLeft; //treshold coefficient LEFT
    private double intensityTresholdRight; //treshold coefficient RIGHT
   
    private double lorentzZeroShift; //LORENTZIAN zero level [p(x)=0]
    private int triangleShift = 5 ; //correction
    
    /**
     * Loads intensity values from txt file into variable intensities.
     */
    
    
     // --- DATA LOADING ---
    public void loadValues(String name) {
        File file = new File(name);
        BufferedReader reader = null;

        try {
            reader = new BufferedReader(new FileReader(file));
            String text = null;

            int row = 0;
            while ((text = reader.readLine()) != null) {
                intensities.add(new ArrayList<>());
                String[] split = text.split("\\s");
                int column = 0;
                for (String s : split) {
                   intensities.get(row).add(new Pixel(row, column, Double.parseDouble(s))); 
                   column++;
                }
                row++;

            }
        } catch (FileNotFoundException e) {
            e.printStackTrace();
        } catch (IOException e) {
            e.printStackTrace();
        } finally {
            try {
                if (reader != null) {
                    reader.close();
                }
            } catch (IOException e) {
            }
        }
        }
    
    
    // --- TIME VARIABLES ---
    public void getSurface(double t) {
        surfaceY =  (int) (-33.314430744 * Math.pow(Math.E, -t / 8.177549447) - 33.427152961 * 
                Math.pow(Math.E, -t / 8.1775548321) + 955.065412107 + 0.5);
    }
    
    public void getBottom(double t) {
        bottomY = (int) (688.3489076 * Math.pow(Math.E, -t / 1.3130046) + 470.1000535 * Math.pow(Math.E, 
                -t / 7.3354304) + 994.1595543 + 0.5); 
    }
    
    public void getSearchLimit(double t) {
        if (t <= 2){
            searchLimit  = 0.6;
        } else if (t <= 15) { 
            searchLimit =  -0.04 * t + 0.8;
        } else if (t <= 20){
            searchLimit = 0.15;
        } else {
            searchLimit = 0.1;
        }
    }
    
    public void getIntensityTresholdLeft(double t) {
        if (t <= 96){
            intensityTresholdLeft  = 1.554716228331147 + 0.8417198413121476*Math.exp(-t / 6.982013495359394) - 0.3 ;
        } else {
            intensityTresholdLeft = 80 ;
        }
    }
    
    public void getIntensityTresholdRight(double t) {
        if (t <= 96){
            intensityTresholdRight  = 1.241399509293841 + 1.005020222773979*Math.exp(-t / 8.926683282760306) - 0.2;
        } else {
            intensityTresholdRight = 90 ;
        }
    }
    
    public void getLorentzZeroShift(double t) {
        if (t <= 7){
        lorentzZeroShift = 41.627914340572 - 74.690983466583 * Math.exp(-t*25/112.20495444165); 
        } else if (t <= 12) {
        lorentzZeroShift = 25;
        } else {
        lorentzZeroShift = 30;
        }
        System.out.println("Zero Shift: " + lorentzZeroShift);
    }
        
    
    // --- POSITION OF VORTEX PEAK (x) ---
    public void findMiddle() {
        double sumInt = 0;
        double sumXInt = 0;
        for (int y = bottomY + 5 ; y > bottomY - 25; y--) {
            for (int x = middleOfContainer - 60; x <= middleOfContainer + 60; x++) {
                  double intensity = intensities.get(y).get(x).value;
                  sumInt += Math.pow(intensity, 4);
                  sumXInt += x * Math.pow(intensity, 4);
            }
        }
        middleOfWhirl = (int) (sumXInt / sumInt + 0.5);
        //System.out.println("peak centroid: " + middleOfWhirl);
    }
    
    /**
     * Stores two transitions from each row into variable transitions
     */
    
    
    // --- VORTEX DETECTION ---
    public void findWhirl() {
        //average background intensity LEFT/RIGHT
        double sumLeft = 0;
        double sumRight = 0;
        
        for (int y = 1000; y < 1040; y++) {
            for (int x = 140; x <= 260; x++) { 
                sumLeft += intensities.get(y).get(x).value;
            }
        }  
        for (int y = 1460; y < 1508; y++) {
            for (int x = 880; x < 980; x++) {
                sumRight += intensities.get(y).get(x).value;
            }
        }
        
        double averageLeft = sumLeft / 4800;
        double averageRight = sumRight / 4800;
        //System.out.println("Left BGR: " + averageLeft);
        //System.out.println("Right BGR: " + averageRight);
        
        bottomY = bottomY + triangleShift; //correction, shifts searching triangle [down]
        
        //left side of triangle - looking for vortex boundary
        for (int y = bottomY ; y >= bottomY - (double)(bottomY - surfaceY)*searchLimit; y--) {
            double param = (double)(middleOfWhirl - leftBoundary) / (bottomY - surfaceY);
            int xBoundary = (int) (param * y - param * surfaceY + leftBoundary + 0.5);
            
            if (y <= bottomY) {
                for (int x = xBoundary; x <= middleOfWhirl; x++) {
                    if (intensities.get(y).get(x).value > averageLeft * intensityTresholdLeft) {
                        transitions.add(new Pixel(x, y, 255));
                        break;
                    }
                }
            }
            
        }
        
        //right side of triangle - looking for vortex boundary
        for (int y = bottomY ; y >= bottomY - (double)(bottomY - surfaceY)*searchLimit; y--) {
            double param = (double)(rightBoundary - middleOfWhirl) / (bottomY - surfaceY);
            int xBoundary = (int) (-param * y + param * surfaceY + rightBoundary + 0.5);
            
            if (y <= bottomY) {
                for (int x = xBoundary; x >= middleOfWhirl; x--) {
                    if (intensities.get(y).get(x).value > averageRight * intensityTresholdRight) {
                        transitions.add(new Pixel(x, y, 255));
                        break;
                    }
                }
            }
            
        }
        
        
    }
        
    
    // --- NONLINEAR REGRESSION: GET PARAMETERS---
    public void fitLorentzian() {
        double[] xdata = new double[transitions.size()];//+ 2];
        double[] ydata = new double[transitions.size()];//+ 2];
        //double[] xerror = new double[transitions.size()];// + 2]; [not used]
        //double[] yerror = new double[transitions.size()];// + 2]; [not used]
        
        for (int i = 0; i < transitions.size(); i++) {
            xdata[i] = transitions.get(i).x; //getting POSITION of"good pixels"
            ydata[i] = transitions.get(i).y - surfaceY - lorentzZeroShift; // !!! SHIFT FOR FIT !!!
        }
        
        // @copyright Dr. Michael Thomas Flanagan
        Regression reg = new Regression(xdata, ydata);// WITH ERRORS: new Regression(xdata, ydata, xerror, yerror);
        reg.lorentzian();
        //reg.gaussian(); //GAUSSIAN

        //reg.print("lorentzfit.txt"); //makes txt file with all the information about regression
        
        double[] estimates = reg.getBestEstimates();
        mean = estimates[0];
        gamma = estimates[1];
        //sigma = estimates[1]; //if GAUSSIAN
        yscale = estimates[2];
        /*for (double param : estimates) {
            System.out.println(param); //prints the estimates from the fit to the output consolee
        }
        */
        System.out.println("Cauchy / Lorentz Function");
        System.out.println("mean: " + mean);
        System.out.println("gamma: " + gamma);
        System.out.println("A0: " + yscale);
    }
    
    
     // --- NONLINEAR REGRESSION: GET CAUCHY CURVE ---
    public void makeEstimation(int number) {
        for (int x = leftBoundary; x <= rightBoundary; x++) {
            transitions.add(new Pixel(x, calculateLorentzFit(x), 255));
        }
        Writer writer = null;

    try {
        writer = new BufferedWriter(new OutputStreamWriter(
        new FileOutputStream("outputPixels".concat(Integer.toString(number)).concat(".txt")), "utf-8")); //creates txt file: "good pixels" + lorentzian fit curve
                                                         
        for (Pixel p : transitions) {
            writer.write(p.x + " " + p.y + "\n");
        }
} catch (IOException ex) {
    // Report
} finally {
   try {writer.close();
   } catch (IOException ex) {/*ignore*/}
}

    }
    
    
    // --- NONLINEAR REGRESSION: CAUCHY FUNCTION FROM PARAMETERS ---
    public int calculateLorentzFit(int x) {
        //double exponent = Math.pow((x-mean)/sigma, 2); //if GAUSSIAN
        //double y = (yscale / (sigma * Math.sqrt(2 * Math.PI)))*Math.exp(-0.5*exponent); //if GAUSSIAN
        double y = (yscale / Math.PI) * ((gamma / 2) / (Math.pow(x - mean, 2) + Math.pow(gamma / 2, 2))); //LORENTZIAN
        //System.out.println("lorentz " + x + " " + y);
        return (int) (y + surfaceY + lorentzZeroShift); // SHIFT FOR FIT (back)
    }
    
    
    // --- EXECUTE PROCESS ---
    public void processImages() {
        double t = 1; // !!! ATTENZIONE !!! first image should be in normal world t=0, here it is t=1 
        for (int i = 1; i <= numberOfImages; i++) { // i...number of image
            System.out.println("Image: " + i);
            System.out.println("Time [s] (+1 s): " + t);
            loadValues("../fast_I_SF10_txt/" + Integer.toString(i).concat(".txt"));
            getSurface(t); //actual position of water surface
            getBottom(t); //actual position of vortex peak (y)
            findMiddle(); //actual position of vortex peak (x)
            getSearchLimit(t); //actual height of search limit
            getIntensityTresholdLeft(t); //actual treshold coefficient LEFT
            getIntensityTresholdRight(t); //actual treshold coefficient RIGHT
            getLorentzZeroShift(t); //actual lorentzian zero level
            t+=0.4; // t increased by the time between two images [according to SAMPLING FREQUENCY]
            findWhirl(); //vortex boundary detection
            fitLorentzian(); //nonlinear regression: parameters of curve
            makeEstimation(i); //nonlinear regression: curve from parameters 
            intensities.clear();
            transitions.clear();
        }
   
    }

    
    public static void main(String[] args)
    {    
        KuboVir whirl = new KuboVir();
        whirl.processImages();
        


           
            //System.out.print("Riadok:" + y + " ");
            
    }     
    
}