// Last edited on 2013-01-05 01:21:35 by stolfilocal
package thog;

import java.awt.BasicStroke;
import java.awt.Color;
import java.awt.Graphics2D;
import java.awt.Polygon;
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;

import javax.imageio.ImageIO;

import JKernelMachines.Classifier;

import utils.image_functions;
import utils.file_functions;
import utils.FloatImage;

/* ----------------------------------------------------------------------
An instance of this class is an extractor of T-HOG image descriptors with 
specified parameters. */
public class thog_desc {
	
    /* Private parameters: */
    public final int new_height;
    public final int number_of_cells_x;
    public final int number_of_cells_y;
    public final int bins_per_cell;
    public final boolean image_normalization;
    public final String image_normalization_weight;
    public final double image_normalization_weight_radius;
    public final int histogram_normalization_metric;
    public final String weight_function;
    public final boolean deformable_weights;
    public final boolean oriented_gradients;
    
    /* === CONSTRUCTORS ===================================================== */

    public thog_desc (String params_file_name) {
        ArrayList<String> thog_params = file_functions.get_list(params_file_name);
        new_height = Integer.parseInt(thog_params.get(0));
        number_of_cells_x = Integer.parseInt(thog_params.get(1));
        number_of_cells_y = Integer.parseInt(thog_params.get(2));
        bins_per_cell = Integer.parseInt(thog_params.get(3));
        image_normalization = Boolean.parseBoolean(thog_params.get(4));
        image_normalization_weight = thog_params.get(5);
        image_normalization_weight_radius = Double.parseDouble(thog_params.get(6));
        String gradient_option = thog_params.get(7); /* !!! ignored for now !!! */
        String norm_metric = thog_params.get(8);
        if (norm_metric.equals("L1")) {
            histogram_normalization_metric = 1;
        } else if (norm_metric.equals("L2")) {
            histogram_normalization_metric = 2;
        } else {
            assert false : "invalid normalization";
            histogram_normalization_metric = 0;
        }
        weight_function = thog_params.get(9);
        deformable_weights = Boolean.parseBoolean(thog_params.get(10));
        /* !!! There should be a parameter to choose oriented vs unoriented gradients. !!! */ 
        oriented_gradients = true;
    }

    /* Extracts a sub-mage with extra margin to ensure that the top and bottom parts of the roman-text are enclosed. */
    public FloatImage get_thog_sub_image (FloatImage img, int x, int y, int w, int h) {
        int margin = (int)Math.floor(0.15*h + 0.5);

        int xmin = x;
        int xmax = x + w - 1;
        int ymin = y - margin;
        int ymax = y + h - 1 + margin;
        
        FloatImage res = img.get_sub_image (xmin, ymin, xmax, ymax);
        return res;
    }	
	
    /* Text/non-text thog classification score. */
    public double score (FloatImage img, Classifier cls, boolean debug, String debug_prefix) {
        double[] descriptor = get_descriptor(img, debug, debug_prefix);
        return cls.valueOf (descriptor);
    }

    /* Extracts the T-HOG descriptor from a given image, after resizing it to a standard height. */
    public double[] get_descriptor (FloatImage img, boolean debug, String debug_prefix) {

        int number_of_cells = number_of_cells_x * number_of_cells_y; 
        int number_of_bins = number_of_cells * bins_per_cell; /*HOG bins*/

        /* Rescale the image to standard height {new_height}. */
        /* If {new_height} is negative the region is not resized. */
        FloatImage gry = null;	 	
        if (new_height > 0) {
            gry = img.scale_to_height(new_height, number_of_cells_x); 
            assert(gry.ny == new_height);
        }
        else {
            gry = img.copy ();
        }
        if (debug){
            /* Writing the extracted and scaled image: */  
            gry.write_as_png (debug_prefix + "gry", (float)0.0, (float)1.0, FloatImage.DEFAULT_GAMMA, true);
        }	

        /*number of pixels of the resized image.*/
        int n = gry.nx * gry.ny;

        /* Margin to exclude from gradient computation: */
        int safe_margin = 1;
            
        /* Computing the approximate baselines {baselines[0..gry.nx-1][0..1]} to each column. */
        double[][] baselines = null;

        if (deformable_weights) {
            /*Computing the variable text-baselines .*/
            int tmp_rwt = choose_gaussian_weight_size (1.0 * gry.ny);
            double[] dist_weight = new double[2 * tmp_rwt + 1];
            compute_gaussian_weights (dist_weight, tmp_rwt);
            baselines = text_baselines (dist_weight, tmp_rwt, gry);
        }
        else {
            /*without variable text-baselines weights*/
            baselines = new double[gry.nx][2];
            for (int x = 0; x < gry.nx; x++) {
                baselines[x][0] = safe_margin; /*ymin*/ 
                baselines[x][1] = gry.ny - safe_margin - 1; /*ymax*/
            }	
        }
        if (debug) {
            /* Write an image showing the estimated baselines: */			
            FloatImage tmp = gry.copy();
            for (int x = safe_margin; x < (gry.nx - safe_margin); x++) {
                /* Get the pixels to plot the baselines extremes.*/
                int ymin = (int)Math.floor(baselines[x][0]);
                int ymax = (int)Math.floor(baselines[x][1]);

                /*writing in black the baselines extremes.*/
                if ((ymin > 0) && (ymin < gry.ny)) {
                    tmp.set_sample(0, x, ymin, -2.0f);
                }
                if ((ymax > 0) && (ymax < gry.ny)) {
                    tmp.set_sample(0, x, ymax, -2.0f);
                }
            }
            tmp.write_as_png (debug_prefix + "baselines", -2.0f, 2.0f, FloatImage.DEFAULT_GAMMA, true);
        }
			
        /* Get images with gradient norm and angle: */
        FloatImage dnorm = new FloatImage(1, gry.nx, gry.ny);
        FloatImage dtheta = new FloatImage(1, gry.nx, gry.ny);
        FloatImage dgrad = ( debug ? new FloatImage(2, gry.nx, gry.ny) : null );
        double noise = 0.03; 
        get_mag_theta (gry, dnorm, dtheta, dgrad, safe_margin, noise);
        if (debug) {
            /* Writing image gradient norm and directions*/			
            dnorm.write_as_png (debug_prefix + "dnorm", 0.0f, 0.5f, FloatImage.DEFAULT_GAMMA, true);
            dtheta.write_as_png (debug_prefix + "dnorm", 0.0f, 0.5f, FloatImage.DEFAULT_GAMMA, true);
            dgrad.write_as_png (debug_prefix + "dgrad", 0.0f, 0.5f, FloatImage.DEFAULT_GAMMA, true);
        }
            
        /*Creating a matrix to hold the cells weights.*/
        FloatImage[][] cell_weight_image = null;
        if (debug) {
            cell_weight_image = new FloatImage[number_of_cells_x][number_of_cells_y];
            for (int cx = 0; cx < number_of_cells_x; cx++) {
               for (int cy = 0; cy < number_of_cells_y; cy++) {
                   cell_weight_image[cx][cy] = new FloatImage(0, gry.nx, gry.ny);
                }
            }
        }

        /* Nominal left and right margins for cell weights: */
        double xleft = safe_margin;
        double xright = gry.nx - safe_margin;

        /*Creating a matrix to hold the cells histogram.*/
        double[] hist = new double[number_of_bins];
        for (int i = 0; i < number_of_bins; i++) {
            hist[i] = 0;
        }
        if ( (gry.nx <= 2) || (gry.ny <= 2) ) {
            return hist;
        }

        /* Accumulate the gradients into the histogram: */
        for (int x = safe_margin; x < (gry.nx - safe_margin); x++) {

            /* Y coordinates of baselines for column {x}.*/
            double ymin = baselines[x][0]; /* Top */
            double ymax = baselines[x][1]; /* Bottom */

            for (int y = safe_margin; y < (gry.ny - safe_margin); y++) {

                /* Get the bin indices where this pixel contributes to: */
                float theta = dtheta.get_sample(0,x,y);
                float norm = dnorm.get_sample(0,x,y);
                int[] bin = new int[2];
                double[] factor = new double[2];
                get_bin_pos (bins_per_cell, theta, bin, factor);
                
                /* Computing the cells histogram and fuzzy weights*/
                for (int cx = 0; cx < number_of_cells_x; cx++) {
                    double wtx = cell_weight (weight_function, number_of_cells_x, cx, x, xleft, xright);
                    for (int cy = 0; cy < number_of_cells_y; cy++) {
                        double wty = cell_weight (weight_function, number_of_cells_y, cy, y, ymin, ymax);
                        /* Index of bin 0 of cell {cx,cy}: */
                        int c_pos = (cy * number_of_cells_x + cx)* bins_per_cell;
                        hist[c_pos + bin[0]] += (norm * wtx * wty) * factor[0];
                        hist[c_pos + bin[1]] += (norm * wtx * wty) * factor[1];
                        if (debug) { 
                            cell_weight_image[cx][cy].set_sample(0, x, y, (float)(wtx * wty));
                        }
                    }
		}
            }
        }

        if (debug) {
            /* writing cell weight functions*/ 
            for (int cx = 0; cx < number_of_cells_x; cx++) {
                for (int cy = 0; cy < number_of_cells_y; cy++) {
                    String cell_tag = "cwt_" + String.format("%02d",cx) + "_" + String.format("%02d",cy);
                    cell_weight_image[cx][cy].write_as_png (debug_prefix + cell_tag, 0.0f, 1.0f, FloatImage.DEFAULT_GAMMA, true);
                }
            }

        }

        /* Normalize the whole histogram to unit L1 or L2 norm: */
        double sum = 0.0;
        for (int cy = 0; cy < number_of_cells_y; cy++) { 
            for (int cx = 0; cx < number_of_cells_x; cx++) {
                /* Index of bin 0 of cell {cx,cy}: */
                int c_pos = (cy * number_of_cells_x + cx) * bins_per_cell;
                for (int bin = 0; bin < bins_per_cell; bin++) {
                    double v = hist[c_pos + bin];
                    if (histogram_normalization_metric == 1) {
                        sum += v; 
                    }
                    else if (histogram_normalization_metric == 2) {
                        sum += v * v;
                    }
                    else {
                        assert false : "invalid normalization";
                    }
                }
            }
        }
        double hist_norm;
        if (histogram_normalization_metric == 1) {
            hist_norm = sum + bins_per_cell*number_of_cells_x*number_of_cells_y;
        }
        else if (histogram_normalization_metric == 2) {
            hist_norm = Math.sqrt(sum + bins_per_cell*number_of_cells_x*number_of_cells_y);
        } else {
            assert false : "invalid normalization";
            hist_norm = 1;
        }
        for (int cy = 0; cy < number_of_cells_y; cy++) { 
            for (int cx = 0; cx < number_of_cells_x; cx++) {
                /* Index of bin 0 of cell {cx,cy}: */
                int c_pos = (cy * number_of_cells_x + cx) * bins_per_cell;
                for (int bin = 0; bin < bins_per_cell; bin++) {
                    int bin_pos = c_pos * bins_per_cell + bin;
                    double v = hist [c_pos + bin];
                    hist [c_pos + bin] = (float)(v/hist_norm);
                }
            }
        }
			
        if (debug) {
            /* Writing hog descriptors as circular images: */
            /* Find range of values over all histograms: */
            double hmax = 0;
            for (int k = 0; k < number_of_bins; k++) { hmax = Math.max(hmax, hist[k]); }
            
            /* Plot each histogram: */
            for (int cy = 0; cy < number_of_cells_y; cy++) {
                for (int cx = 0; cx < number_of_cells_x; cx++) {
                    /* Extract the histogram of cell {cx,cy}: */
                    int c_pos = (cy * number_of_cells_x + cx) * bins_per_cell;
                    double[] histogram = Arrays.copyOfRange(hist, c_pos, c_pos + bins_per_cell - 1);
                    String file_name = debug_prefix + "hog_directions_" + cy + "_" + cx + ".png";
                    print_histogram (file_name, histogram, bins_per_cell, hmax);
		}
            }
        }
        
        return hist;
    } 	

    public void get_bin_pos (int bins_per_cell, double dtheta, int[] bin, double[] factor) {
        /*Computing the bin according the gradient direction.*/

        double period = (oriented_gradients ? 2 * Math.PI : Math.PI);

        double a = bins_per_cell*(dtheta/period+1);
        int ia = (int)Math.floor(a);
        double frac = a - (double)ia;

        bin[0] = (ia + bins_per_cell) % bins_per_cell;
        bin[1] = (ia + 1) % bins_per_cell;
        factor[0] = (1-frac);
        factor[1] = frac;

        assert ( (bin[0] >= 0) && (bin[0] < bins_per_cell));
        assert ( (bin[1] >= 0) && (bin[1] < bins_per_cell));		
    }

    /* ----------------------------------------------------------------------
    Computes the norm and angle of the gradient of image {gry}, stores them in {dnorm,dtheta}.  If {dgrad} is not null,
    also saves the gradient components as channels 0 and 1 of {dgrad}. Skips {safe_margin} pixels
    around the border. Kills the gradient norm if it is less than {noise*sqrt(2)}.
    The angle ranges in {[0 _ PI]} or {[0 _ 2*PI]} depending on {oriented_gradients}. */ 
    public void get_mag_theta (FloatImage gry, FloatImage dnorm, FloatImage dtheta, FloatImage dgrad, int safe_margin, double noise) {

        double eps = Math.sqrt(2.0)*noise; /*Assumed deviation of noise in gradient*/
        double eps2 = eps * eps;
        for (int y = safe_margin; y < (gry.ny - safe_margin); y++) {
            for (int x = safe_margin; x < (gry.nx - safe_margin); x++) {
                /*Computing image gradients*/
                double[] grad = gry.gradient_simple(0, x, y);
                if (dgrad != null) {
                    dgrad.set_sample(0, x, y, (float)grad[0]);
                    dgrad.set_sample(1, x, y, (float)grad[1]);
                }
                /*Computing the gradient norm but return zero if too small*/
                double d2 = grad[0]*grad[0] + grad[1]*grad[1];
                float norm;
                if (d2 <= eps2) {
                    norm = 0.0f;
                } else {
                    norm = (float)Math.sqrt(d2 - eps2);
                }
		dnorm.set_sample(0, x, y, norm);

                /* Computing the gradient direction.*/
                float theta = (float)Math.atan2(grad[1], grad[0]);
                if (theta < 0) { theta += (oriented_gradients ? 2*Math.PI : Math.PI); }
            }
	}
    }

    /* Computes the cell weight factor for one axis {z} (x or y). Given: number of cells 
     * {ncz}, cell index {cz}, pixel index {z}, nominal range {zmin,zmax}, all along that axis. 
     * The option selects the weight type. Adjusting the interval by {r_bot, r_top}.*/
    public double cell_weight (String weight_function, int ncz, int cz, int z, double zmin, double zmax) {

        if (ncz == 1) { return 1; }

        double z_star = (z - zmin)/(zmax - zmin);

        if (weight_function.compareTo("Step") == 0) {
            return StepFunc (ncz, cz, z_star);
        } else if (weight_function.compareTo("Bernstein") == 0){
            return Bernstein (ncz - 1, cz, z_star);
        } else if (weight_function.compareTo("Core") == 0) {
            return EdgeCore (ncz, cz, z_star);
        } else if (weight_function.compareTo("Exp") == 0){
            return Exp (ncz, cz, z_star);
        } else {
            assert false : "invalid weight function";
            return 1.0;
	}
    }

    /*Divides the interval [0-1] into {n} equal parts and returns 1.0 if {z} 
     *is in part number {k} (0..n-1), 0 otherwise. */
    public double StepFunc (int n, int k, double z) {
        assert((k >= 0) && (k < n));
        if (z <= 0) { return (k == 0 ? 1 : 0); }
        else if (z >= 1) { return (k == (n - 1) ? 1 : 0); }
        else {
            return ( (k <= z*n) && (z*n < k+1) ? 1 : 0);
        }
    }

    /*Computes the Bernstein polynomial of degree {n} and index {k} for the argument {z}.*/
    public double Bernstein (int n, int k, double z) {
        assert((k >= 0) && (k <= n));
        if (z <= 0) { return (k == 0 ? 1 : 0); }
        else if (z >= 1) { return (k == n ? 1 : 0); }
        else {
            double zmax = ((double)k)/((double)n);
            return BernsteinPoly(n,k,z)/BernsteinPoly(n,k,zmax);
        }
    }

    public double BernsteinPoly (int n, int k, double z) {
        assert((k >= 0) && (k <= n));
        double res = 1.0;
        for (int i = 0; i < k; i++) {
            res = (res * (n - i))/(i+1)*z;
        }
        return res*Math.pow(1-z,n-k);
    }

    /*An edge-core weight function. If {n == 1} returns 1, if (n == 2) returns 
     *weight 1.0 near the edges, or 1.0 in the core region depending on {k}*/
    public double EdgeCore (int n, int k, double z) {
        assert(n == 2);
        assert((k >= 0) && (k < n));
        if ( (z <= 0) || (z >= 1) ) { return (k == 0 ? 1 : 0); }
        else {
            double v = 4 * z * (1 - z);
            v = v*v;
            return (k==0? 1 - v: v);
        }
    }

    /*An exponential weight function.*/
    double Exp (int n, int k, double z) {

        double mu = 0.01; 

        double sigma = 0.5; 

        assert((k >= 0) && (k < n));

        double avg = -mu + (1 + 2 * mu)/(double)(n - 1)*k;

        double dev = sigma/n;

        return gaussian (z, avg, dev);		
    }
	
    private double gaussian (double z, double mu, double sigma) {
        return Math.exp(- ((z - mu)*(z - mu))/(2 * sigma * sigma));
    }

    /**/
    public int choose_gaussian_weight_size (double dev) {
        double rwt = 2*dev;
        return (int)(Math.ceil(rwt));
    }

    /**/
    public void compute_binomial_weights (double[] weight, int rwt) {
        int nwt = 2*rwt+1;
        weight[0] = 1;
        for (int i = 1; i < nwt; i++) { 
            weight[i] = 0.0; 
            for (int j = i; j >=1; j--) {
                weight[j] = (weight[j] + weight[j-1])/2;
            }
            weight[0] /= 2;
        }
    }

    /**/
    public void compute_gaussian_weights (double[] weight, int rwt) {
        double eps = 0.01;
        double a = -Math.log(eps);
        int nwt = 2*rwt+1;
        
        for (int i = 0; i < nwt; i++) { 
            double z = (i - rwt)/((double)rwt);
            weight[i] = (1-z*z)*Math.exp(-a*z*z);
        }
    }

    /*************************Drawing HOG histograms**************************/

    /*Function used to draw the HOG descriptor with a circular shape, 
     *where each bin is proportional to its gradient norm.*/
    public void print_histogram (String outname, double[] histogram, int bins_per_cell, double hmax) {

        /* HOG descriptor image: */
        int nx = 600;		
        int ny = 600;

        /* Inner circle radius: */
        double in_radius = 90;
        double out_radius = 300;

        /* Diagram center: */
        double c_x = nx/2;		
        double c_y = ny/2;

	BufferedImage temp = new BufferedImage(nx, ny, BufferedImage.TYPE_INT_RGB);
        /*Defining the graphics settings.*/
        float[] hsbvals = new float[3];		
        Graphics2D g = temp.createGraphics();		
        g.setStroke(new BasicStroke(1.0f));			
        Color.RGBtoHSB(0, 0, 0, hsbvals);		
        g.setColor(Color.getHSBColor(hsbvals[0], hsbvals[1], hsbvals[2]));

        /*Background painting*/		
        int white = 0xFFFFFF;
        for (int y = 0; y < ny; y++) {
            for (int x = 0; x < ny; x++) {
                temp.setRGB(x, y, white);
            }
        }

	double period = (oriented_gradients ? 2 * Math.PI : Math.PI);

        for (int bin = 0; bin < bins_per_cell; bin++) {
            /* Nominal angle range for this bin: */
            double theta_a = (bin+0.4)*period/bins_per_cell;
            double theta_b = (bin-0.4)*period/bins_per_cell;
            /* Mass in this bin: */
            double hb = histogram[bin]/hmax;

            Polygon polygon = null;
            polygon = draw_Polygon (c_x, c_y, in_radius, out_radius, theta_a, theta_b, hb);
            g.fillPolygon(polygon);
            
            if (! oriented_gradients) {
                polygon = draw_Polygon (c_x, c_y, in_radius, out_radius, theta_a + Math.PI, theta_b + Math.PI, hb);
                g.fillPolygon(polygon);
            }
        }
        g.dispose();
        try {
            ImageIO.write(temp, "png", new File(outname));
        }
        catch (Exception e) { 
            System.err.printf("cannot write image\n");
        }
    }

    /* Function to draw each HOG bin as a trapezoid around a circle.*/
    /* Maps {hb} from {[0 _ 1]} to {[in_radius _ out_radius]}. */
    /* Assumes that the Y axis points up.*/
    private Polygon draw_Polygon (double c_x, double c_y, double in_radius, double out_radius, double theta_a, double theta_b, double hb) { 
        double dr = out_radius - in_radius;
        
        /*Taking each point of the rectangle*/
        double x1 = c_x + in_radius * Math.cos(theta_a);
        double y1 = c_y + in_radius * Math.sin(theta_a);

        double x2 = c_x + in_radius + dr*hb * Math.cos(theta_a);
        double y2 = c_y + in_radius + dr*hb * Math.sin(theta_a);

        double x3 = c_x + in_radius * Math.cos(theta_b);
        double y3 = c_y + in_radius * Math.sin(theta_b);

        double x4 = c_x + in_radius + dr*hb * Math.cos(theta_b);
        double y4 = c_y + in_radius + dr*hb * Math.sin(theta_b);

        /*Getting the polygon*/
        Polygon polygon = new Polygon();
        polygon.addPoint((int)x1, (int)y1);
        polygon.addPoint((int)x2, (int)y2);
        polygon.addPoint((int)x4, (int)y4);
        polygon.addPoint((int)x3, (int)y3); 

        return polygon;
    }

    /* Returns an array {baselines[0..nx-1][0..1]} where {baselines[x][0]} is the estimated top baseline
    of column {x}, {baselines[x][1]} is the estimated bottom baseline. */
    public double[][] text_baselines (double[] weight, int rwt, FloatImage gry) {

        double[] s = new double[gry.ny];
        double[][] baselines = new double[gry.nx][2]; 

        for (int x = 0; x < gry.nx; x++) {

            for (int y = 0; y < gry.ny; y++) {	
                double sum_var = 0.0, sum_avg = 0.0, sum_w = 0.0;
                for (int dx = -rwt; dx <= rwt; dx++) {
                    int x1 = x + dx;
                    if ((x1 >= 0) && (x1 < gry.nx)) {
                        double w = weight[rwt+dx];
                        double v = gry.get_sample(0,x,y);
                        sum_avg += w*v;
                        sum_w += w; 
                    }
                }	
                assert(sum_w != 0.0); 

                double AVG = (sum_avg/sum_w);
                
                for (int dx = -rwt; dx <= rwt; dx++) {
                    int x1 = x + dx;
                    if ((x1 >= 0) && (x1 < gry.nx)) {
                        double w = weight[rwt+dx];
                        double v = gry.get_sample(0,x,y);
                        sum_var += w*(v-AVG)*(v-AVG);
                    }
                }				
                s[y] = (sum_var/sum_w);
            }
            double[] tmp = get_baselines_by_percentiles (gry.ny, s);
            baselines[x][0] = tmp[0];
            baselines[x][1] = tmp[1];
        }
        return baselines;
    }

    public double[] get_baselines_by_percentiles (int ny, double[] s) {
        double percentile = 0.2; 
        double[] baselines = new double[2];

        double sum_s = 0.0;
        for (int y = 0; y < ny; y++) { sum_s += s[y]; }

        int top = 0;
        double sum_top = 0.0;
        while ( (top < ny) && (sum_top+s[top] <= percentile*sum_s) ) {
            sum_top += s[top];
            top++;
        }
        if (top > 0) { top--; }

        int bot = ny-1;
        double sum_bot = 0.0;
        while ( (bot >= 0) && (sum_bot+s[bot] <= percentile*sum_s) ) {
            sum_bot += s[bot];
            bot--;
        }
        if (bot < ny-1) { bot++; }

        double y_med = 0.5*(top+bot);
        double y_delta_min = 0.2*ny;
        double y_delta = Math.max(0.5*(bot-top), y_delta_min);

        baselines[0] = Math.max(y_med - y_delta,0);   //ymin, TOP 
        baselines[1] = Math.min(y_med + y_delta,ny-1); //ymax, BOTTOM
        
        return baselines;
    }

    public double[] get_baselines_by_med (int ny, double[] s) {
        double[] baselines = new double[2];
        double sum_s = 0.0;
        double sum_sy = 0.0;

        for (int y = 0; y < ny; y++) {
            sum_sy += s[y] * y;
            sum_s += s[y];
        }
        
        double y_med = sum_sy/sum_s;
        double sum_sdy2 = 0.0;
        for (int y = 0; y < ny; y++) {
            sum_sdy2 += s[y] * ( (y - y_med)*(y - y_med) );
            sum_s += s[y];
        }

        double y_dev = Math.sqrt(sum_sdy2/sum_s);

        double y_delta_min = 0.2*ny;

        double y_delta = Math.max(3*y_dev, y_delta_min);

        baselines[0] = Math.max(y_med - y_delta,0); /*ymin, TOP*/ 
        baselines[1] = Math.min(y_med + y_delta,ny); /*ymax, BOTTOM*/

        return baselines;
    }
}